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=(0,10),skiprows=10,unpack=True)
  33  x1,e,r,gy,nc= np.loadtxt('md.stat',usecols=(0,3,5,12,6),skiprows=0,unpack=True)
  34 else:
  35  x,y= np.loadtxt('md.stat',usecols=(0,6),skiprows=10,unpack=True)
  36  x1,e,gy= np.loadtxt('md.stat',usecols=(0,3,8),skiprows=0,unpack=True)
  37
  38 h,bin=np.histogram(y,bins=50,density=True)
  39
  40 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
  41 plt.xlim(min(bin), max(bin))
  42 plt.ylim(0,max(h))
  43 plt.ylabel('probality')
  44 plt.xlabel('temperature')
  45
  46 center = (bin[:-1] + bin[1:]) / 2
  47 #print bin
  48 #print center
  49 start = (g-2)/2*np.log(t_bath) - g*t_bath/(2*t_bath)
  50 popt, pcov = curve_fit(prob_T, center, h, p0=-start)
  51 xfine = np.linspace(min(bin), max(bin), 100)
  52 #print popt
  53
  54 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
  55 print '%15.10f' % chi_squared
  56
  57 #print xfine
  58 #print prob_T(xfine,popt[0])
  59 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
  60 plt.savefig('temp_hist.png')
  61 #plt.show()
  62
  63 plt.clf()
  64 plt.xlabel('step')
  65 plt.ylabel('potential energy')
  66 plt.plot(x1,e,'.')
  67 plt.savefig('md_ene.png')
  68
  69 plt.clf()
  70 plt.xlabel('step')
  71 plt.ylabel('radius of gyration')
  72 plt.plot(x1,gy,'.')
  73 plt.savefig('md_gyr.png')
  74
  75 if ncolumns==14:
  76  plt.clf()
  77  plt.xlabel('step')
  78  plt.ylabel('RMSD')
  79  plt.plot(x1,r,'.')
  80  plt.savefig('md_rms.png')
  81
  82  plt.clf()
  83  plt.xlabel('step')
  84  plt.ylabel('fraction of native side-chain concacts')
  85  plt.plot(x1,nc,'.')
  86  plt.savefig('md_fracn.png')
  87