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=40,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=40,range=(min(binall),max(binall[hall>4])),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 cmap = plt.get_cmap('hot')
  67 colors = cmap(np.linspace(0, 1, len(replica)*1.4))
  68 for i,c in zip(replica,colors):
  69  yt=x0[r0==i]
  70  xt=(s0+r0*max(s0))[r0==i]
  71  plt.plot(xt,yt,'-',color=c)
  72
  73 plt.savefig('remd_ex.png')
  74
  75 colors = cm.rainbow(np.linspace(0, 1, len(Tremd)))
  76 if ncolumns==14:
  77
  78   plt.clf()
  79   plt.xlabel('rmsd')
  80   plt.ylabel('potential energy')
  81
  82   for T,c in zip(Tremd,colors):
  83     xt=rms[x==T]
  84     yt=y[x==T]
  85     plt.plot(xt,yt,'.',color=c,ms=4)
  86
  87
  88   plt.savefig('remd_ene_rms.png')
  89
  90   plt.clf()
  91   plt.xlabel('step*replica')
  92   plt.ylabel('rmsd')
  93   for i in replica:
  94     yt=rms0[r0==i]
  95     xt=(s0+r0*max(s0))[r0==i]
  96     tt=x0[r0==i]
  97     plt.scatter(xt,yt,c=tt,edgecolors='face',s=0.1,cmap=cm.rainbow,vmin=Tremd[0],vmax=Tremd[-1])
  98   plt.xlim(0,max(s)+max(s)*max(r))
  99   plt.savefig('remd_step_rms.png')
 100
 101
 102 x,y,rms= np.loadtxt('file_wham.thermal',usecols=(0,6,4),unpack=True)
 103
 104 plt.clf()
 105 plt.xlabel('bath temperature [K]')
 106 plt.ylabel('heat capacity')
 107 plt.xlim(Tremd[0]-10, Tremd[-1]+10)
 108 plt.plot(x,y,'-',color=c)
 109 plt.savefig('remd_cv.png')
 110
 111 if ncolumns==14:
 112   plt.clf()
 113   plt.xlabel('bath temperature [K]')
 114   plt.ylabel('average RMSD')
 115   plt.xlim(Tremd[0]-10, Tremd[-1]+10)
 116   plt.plot(x,rms,'-')
 117   plt.savefig('remd_rmsd.png')
 118