6 import matplotlib.pyplot as plt
8 from scipy.optimize import curve_fit
14 aa=np.float128(10**(-gg-2)*a)
16 return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
17 # return aa * ( x**((gg-2)/2) * np.exp( -gg*x/(2*Tr) ) )
19 #x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
20 #x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
21 #x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
23 with open('md.stat','r') as f:
26 ncolumns=len(line.split())
29 x,y= np.loadtxt('md.stat',usecols=(0,10),skiprows=10,unpack=True)
30 x1,e,r,gy,nc= np.loadtxt('md.stat',usecols=(0,3,5,12,6),skiprows=0,unpack=True)
32 x,y= np.loadtxt('md.stat',usecols=(0,6),skiprows=10,unpack=True)
33 x1,e,gy= np.loadtxt('md.stat',usecols=(0,3,8),skiprows=0,unpack=True)
35 h,bin=np.histogram(y,bins=50,density=True)
37 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
38 plt.xlim(min(bin), max(bin))
40 plt.ylabel('probality')
41 plt.xlabel('temperature')
43 center = (bin[:-1] + bin[1:]) / 2
46 popt, pcov = curve_fit(prob_T, center, h)
47 xfine = np.linspace(min(bin), max(bin), 100)
50 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
51 print '%15.10f' % chi_squared
54 #print prob_T(xfine,popt[0])
55 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
56 plt.savefig('temp_hist.png')
61 plt.ylabel('potential energy')
63 plt.savefig('md_ene.png')
67 plt.ylabel('radius of gyration')
69 plt.savefig('md_gyr.png')
76 plt.savefig('md_rms.png')
80 plt.ylabel('fraction of native side-chain concacts')
82 plt.savefig('md_fracn.png')