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  x0,s0,r0,rms0= np.loadtxt('remd_all0.stat',usecols=(11,0,13,5),unpack=True)
  18 else:
  19  x,y,s,r= np.loadtxt('remd_all.stat',usecols=(7,3,0,9),unpack=True)
  20  x0,s0,r0= np.loadtxt('remd_all0.stat',usecols=(7,0,9),unpack=True)
  21
  22 hall,binall=np.histogram(y,bins='auto',density=False)
  23
  24 plt.xlim(min(binall), max(binall[hall>4]))
  25 #plt.ylim(0,max(hall)/4)
  26 plt.ylabel('number of samples')
  27 plt.xlabel('potential energy [kcal/mol]')
  28
  29 #Tremd=[240, 260, 280, 300, 320, 340, 360, 390]
  30
  31 Tremd=map(float,sys.argv[1].split())
  32
  33
  34 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  35 for T,c in zip(Tremd,colors):
  36  yt=y[x==T]
  37  h,bin=np.histogram(yt,bins='auto',density=False)
  38  center = (bin[:-1] + bin[1:]) / 2
  39  plt.plot(center,h,'-',color=c)
  40 # plt.bar(bin[:-1], h, width = bin[2]-bin[1],color=c)
  41
  42 plt.savefig('remd_ene_hist.png')
  43 #plt.show()
  44
  45 plt.clf()
  46 plt.xlabel('bath temperature [K]')
  47 plt.ylabel('potential energy [kcal/mol]')
  48
  49 plt.ylim(min(binall), max(binall[hall>4]))
  50 plt.xlim(Tremd[0]-10, Tremd[-1]+10)
  51 #Tremd=[240, 260, 280, 300, 320, 340, 360, 390]
  52 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  53 for T,c in zip(Tremd,colors):
  54  yt=y[x==T]
  55  xt=x[x==T]
  56  plt.plot(xt,yt,'.',color=c)
  57
  58 plt.savefig('remd_Tene.png')
  59
  60 plt.clf()
  61 plt.ylabel('bath temperature [K]')
  62 plt.xlabel('step*replica')
  63
  64 replica=range(int(sys.argv[2]))
  65 colors = cm.rainbow(np.linspace(0, 1, len(replica)))
  66 for i,c in zip(replica,colors):
  67  yt=x0[r0==i]
  68  xt=(s0+r0*max(s0))[r0==i]
  69  plt.plot(xt,yt,'-',color=c)
  70
  71 plt.savefig('remd_ex.png')
  72
  73 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  74 if ncolumns==14:
  75
  76   plt.clf()
  77   plt.xlabel('rmsd')
  78   plt.ylabel('potential energy')
  79
  80   for T,c in zip(Tremd,colors):
  81     xt=rms[x==T]
  82     yt=y[x==T]
  83     plt.plot(xt,yt,'.',color=c,ms=4)
  84
  85
  86   plt.savefig('remd_ene_rms.png')
  87
  88   plt.clf()
  89   plt.xlabel('step*replica')
  90   plt.ylabel('rmsd')
  91   for i in replica:
  92     yt=rms0[r0==i]
  93     xt=(s0+r0*max(s0))[r0==i]
  94     tt=x0[r0==i]
  95     plt.scatter(xt,yt,c=tt,edgecolors='face',s=0.1,cmap=cm.rainbow,vmin=Tremd[0],vmax=Tremd[-1])
  96   plt.xlim(0,max(s)+max(s)*max(r))
  97   plt.savefig('remd_step_rms.png')
  98
  99
 100 x,y,rms= np.loadtxt('file_wham.thermal',usecols=(0,6,4),unpack=True)
 101
 102 plt.clf()
 103 plt.xlabel('bath temperature [K]')
 104 plt.ylabel('heat capacity')
 105 plt.xlim(Tremd[0]-10, Tremd[-1]+10)
 106 plt.plot(x,y,'-',color=c)
 107 plt.savefig('remd_cv.png')
 108
 109 if ncolumns==14:
 110   plt.clf()
 111   plt.xlabel('bath temperature [K]')
 112   plt.ylabel('average RMSD')
 113   plt.xlim(Tremd[0]-10, Tremd[-1]+10)
 114   plt.plot(x,rms,'-')
 115   plt.savefig('remd_rmsd.png')
 116