files/matplotlib_fit_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 numpy as np
   8 from scipy.optimize import curve_fit
   9 from math import sqrt
  10 import sys
  11
  12 def prob_T(x,a):
  13     gg=np.float128(g)
  14 #    aa=np.float128(10**(-gg-2)*a)
  15     aa=np.float128(a)
  16     Tr=np.float128(t_bath)
  17     return np.exp( aa + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
  18 #    return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
  19 #    return aa * (  x**((gg-2)/2) * np.exp( -gg*x/(2*Tr) ) )
  20
  21 #x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
  22 #x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
  23 #x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
  24 g=int(sys.argv[1])
  25 t_bath=float(sys.argv[2])
  26 with open('md.stat','r') as f:
  27   for line in f:
  28     pass
  29   ncolumns=len(line.split())
  30
  31 if ncolumns==14:
  32  x,y= np.loadtxt('md.stat',usecols=(1,10),skiprows=10,unpack=True)
  33  x1,e,r,gy,nc= np.loadtxt('md.stat',usecols=(1,3,5,12,6),skiprows=0,unpack=True)
  34 else:
  35  x,y= np.loadtxt('md.stat',usecols=(1,6),skiprows=10,unpack=True)
  36  x1,e,gy= np.loadtxt('md.stat',usecols=(1,3,8),skiprows=0,unpack=True)
  37
  38 x=x*0.0489
  39 x1=x1*0.0489
  40
  41 h,bin=np.histogram(y,bins=50,density=True)
  42
  43 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
  44 plt.xlim(min(bin), max(bin))
  45 plt.ylim(0,max(h))
  46 plt.ylabel('probality')
  47 plt.xlabel('temperature')
  48
  49 center = (bin[:-1] + bin[1:]) / 2
  50 #print bin
  51 #print center
  52 start = (g-2)/2*np.log(t_bath) - g*t_bath/(2*t_bath)
  53 popt, pcov = curve_fit(prob_T, center, h, p0=-start)
  54 xfine = np.linspace(min(bin), max(bin), 100)
  55 #print popt
  56
  57 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
  58 print '%15.10f' % chi_squared
  59
  60 #print xfine
  61 #print prob_T(xfine,popt[0])
  62 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
  63 plt.savefig('temp_hist.png')
  64 #plt.show()
  65
  66 plt.clf()
  67 plt.xlabel('time [ps]')
  68 plt.ylabel('potential energy')
  69 plt.plot(x1,e,'.')
  70 plt.savefig('md_ene.png')
  71
  72 plt.clf()
  73 plt.xlabel('time [ps]')
  74 plt.ylabel('radius of gyration')
  75 plt.plot(x1,gy,'.')
  76 plt.savefig('md_gyr.png')
  77
  78 if ncolumns==14:
  79  plt.clf()
  80  plt.xlabel('time [ps]')
  81  plt.ylabel('RMSD')
  82  plt.plot(x1,r,'.')
  83  plt.savefig('md_rms.png')
  84
  85  plt.clf()
  86  plt.xlabel('time [ps]')
  87  plt.ylabel('fraction of native side-chain concacts')
  88  plt.plot(x1,nc,'.')
  89  plt.savefig('md_fracn.png')
  90