6 import matplotlib.pyplot as plt
8 from scipy.optimize import curve_fit
14 # aa=np.float128(10**(-gg-2)*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) ) )
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)
25 t_bath=float(sys.argv[2])
26 with open('md.stat','r') as f:
29 ncolumns=len(line.split())
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)
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)
38 h,bin=np.histogram(y,bins=50,density=True)
40 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
41 plt.xlim(min(bin), max(bin))
43 plt.ylabel('probality')
44 plt.xlabel('temperature')
46 center = (bin[:-1] + bin[1:]) / 2
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)
54 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
55 print '%15.10f' % chi_squared
58 #print prob_T(xfine,popt[0])
59 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
60 plt.savefig('temp_hist.png')
65 plt.ylabel('potential energy')
67 plt.savefig('md_ene.png')
71 plt.ylabel('radius of gyration')
73 plt.savefig('md_gyr.png')
80 plt.savefig('md_rms.png')
84 plt.ylabel('fraction of native side-chain concacts')
86 plt.savefig('md_fracn.png')