X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=files%2Fmatplotlib_fit_hist.py;h=b31025e51bdfc67879fdf30e65ebc743177b4fba;hb=d03bf58e6dc3d690486f8eb7a8f18ad34938c71f;hp=a932a04d23d1e2309cb696d5b596cb0d6f6a44dc;hpb=591ac48194b207c02ffd7cda59c1c9709a114498;p=django_unres.git

diff --git a/files/matplotlib_fit_hist.py b/files/matplotlib_fit_hist.py
index a932a04..b31025e 100755
--- a/files/matplotlib_fit_hist.py
+++ b/files/matplotlib_fit_hist.py
@@ -11,26 +11,32 @@ import sys
 
 def prob_T(x,a):
     gg=np.float128(g)
-    aa=np.float128(10**(-gg-2)*a)
-    Tr=np.float128(300.)
-    return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
+#    aa=np.float128(10**(-gg-2)*a)
+    aa=np.float128(a)
+    Tr=np.float128(t_bath)
+    return np.exp( aa + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
+#    return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
 #    return aa * (  x**((gg-2)/2) * np.exp( -gg*x/(2*Tr) ) )
 
 #x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
 #x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
 #x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
 g=int(sys.argv[1])
-with open('file_GB000.stat','r') as f:
+t_bath=float(sys.argv[2])
+with open('md.stat','r') as f:
   for line in f:
     pass
   ncolumns=len(line.split())
 
 if ncolumns==14:  
- x,y= np.loadtxt('file_GB000.stat',usecols=(0,10),skiprows=10,unpack=True)
- x1,e,r= np.loadtxt('file_GB000.stat',usecols=(0,3,5),skiprows=0,unpack=True)
+ x,y= np.loadtxt('md.stat',usecols=(1,10),skiprows=10,unpack=True)
+ x1,e,r,gy,nc= np.loadtxt('md.stat',usecols=(1,3,5,12,6),skiprows=0,unpack=True)
 else:
- x,y= np.loadtxt('file_GB000.stat',usecols=(0,6),skiprows=10,unpack=True)
- x1,e= np.loadtxt('file_GB000.stat',usecols=(0,3),skiprows=0,unpack=True)
+ x,y= np.loadtxt('md.stat',usecols=(1,6),skiprows=10,unpack=True)
+ x1,e,gy= np.loadtxt('md.stat',usecols=(1,3,8),skiprows=0,unpack=True)
+
+x=x*0.0489
+x1=x1*0.0489
  
 h,bin=np.histogram(y,bins=50,density=True)
 
@@ -43,7 +49,8 @@ plt.xlabel('temperature')
 center = (bin[:-1] + bin[1:]) / 2
 #print bin
 #print center
-popt, pcov = curve_fit(prob_T, center, h)
+start = (g-2)/2*np.log(t_bath) - g*t_bath/(2*t_bath)
+popt, pcov = curve_fit(prob_T, center, h, p0=-start)
 xfine = np.linspace(min(bin), max(bin), 100)
 #print popt
 
@@ -57,15 +64,27 @@ plt.savefig('temp_hist.png')
 #plt.show()   
 
 plt.clf()
-plt.xlabel('step')
+plt.xlabel('time [ps]')
 plt.ylabel('potential energy')
 plt.plot(x1,e,'.')
 plt.savefig('md_ene.png')
 
+plt.clf()
+plt.xlabel('time [ps]')
+plt.ylabel('radius of gyration')
+plt.plot(x1,gy,'.')
+plt.savefig('md_gyr.png')
+
 if ncolumns==14:
  plt.clf()
- plt.xlabel('step')
+ plt.xlabel('time [ps]')
  plt.ylabel('RMSD')
  plt.plot(x1,r,'.')
  plt.savefig('md_rms.png')
- 
\ No newline at end of file
+
+ plt.clf()
+ plt.xlabel('time [ps]')
+ plt.ylabel('fraction of native side-chain concacts')
+ plt.plot(x1,nc,'.')
+ plt.savefig('md_fracn.png')
+