files/matplotlib_hist.py

   1 #! /usr/bin/env python
   2
   3 import matplotlib
   4 #matplotlib.use('GTK')
   5 matplotlib.use('Agg')
   6 import matplotlib.pyplot as plt
   7 import matplotlib.cm as cm
   8 import numpy as np
   9 import sys
  10
  11 with open('remd_all.stat','r') as f:
  12   line=f.readline()
  13   ncolumns=len(line.split())
  14
  15 if ncolumns==14:
  16  x,y,s,r,ek,rms= np.loadtxt('remd_all.stat',usecols=(11,3,0,13,2,5),unpack=True)
  17 else:
  18  x,y,s,r= np.loadtxt('remd_all.stat',usecols=(7,3,0,9),unpack=True)
  19
  20 hall,binall=np.histogram(y,bins=40,density=False)
  21
  22 plt.xlim(min(binall), max(binall[hall>4]))
  23 plt.ylim(0,max(hall)/4)
  24 plt.ylabel('number of samples')
  25 plt.xlabel('potential energy [kcal/mol]')
  26
  27 #Tremd=[240, 260, 280, 300, 320, 340, 360, 390]
  28
  29 Tremd=map(float,sys.argv[1].split())
  30
  31
  32 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  33 for T,c in zip(Tremd,colors):
  34  yt=y[x==T]
  35  h,bin=np.histogram(yt,bins=40,density=False)
  36  center = (bin[:-1] + bin[1:]) / 2
  37  plt.plot(center,h,'-',color=c)
  38 # plt.bar(bin[:-1], h, width = bin[2]-bin[1],color=c)
  39
  40 plt.savefig('remd_ene_hist.png')
  41 #plt.show()
  42
  43 plt.clf()
  44 plt.xlabel('bath temperature [K]')
  45 plt.ylabel('potential energy [kcal/mol]')
  46
  47 plt.ylim(min(binall), max(binall[hall>4]))
  48 plt.xlim(Tremd[0]-10, Tremd[-1]+10)
  49 #Tremd=[240, 260, 280, 300, 320, 340, 360, 390]
  50 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  51 for T,c in zip(Tremd,colors):
  52  yt=y[x==T]
  53  xt=x[x==T]
  54  plt.plot(xt,yt,'.',color=c)
  55
  56 plt.savefig('remd_Tene.png')
  57
  58 plt.clf()
  59 plt.ylabel('bath temperature [K]')
  60 plt.xlabel('step*replica')
  61
  62 replica=range(int(sys.argv[2]))
  63 colors = cm.rainbow(np.linspace(0, 1, len(replica)))
  64 for i,c in zip(replica,colors):
  65  yt=x[r==i]
  66  xt=(s+r*max(s))[r==i]
  67  plt.plot(xt,yt,'-',color=c)
  68
  69 plt.savefig('remd_ex.png')
  70
  71 if ncolumns==14:
  72
  73   plt.clf()
  74   plt.xlabel('rmsd')
  75   plt.ylabel('potential energy')
  76
  77   for T,c in zip(Tremd,colors):
  78     xt=rms[x==T]
  79     yt=y[x==T]
  80     plt.plot(xt,yt,'.',color=c,ms=4)
  81
  82
  83   plt.savefig('remd_ene_rms.png')
  84
  85
  86
  87 x,y,rms= np.loadtxt('file_wham.thermal',usecols=(0,6,4),unpack=True)
  88
  89 plt.clf()
  90 plt.xlabel('bath temperature [K]')
  91 plt.ylabel('heat capacity')
  92 plt.xlim(Tremd[1]-10, Tremd[-1]+10)
  93 plt.plot(x,y,'-',color=c)
  94 plt.savefig('remd_cv.png')
  95
  96 if ncolumns==14:
  97   plt.clf()
  98   plt.xlabel('bath temperature [K]')
  99   plt.ylabel('average RMSD')
 100   plt.xlim(Tremd[1]-10, Tremd[-1]+10)
 101   plt.plot(x,rms,'-')
 102   plt.savefig('remd_rmsd.png')