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