ctest multichain protein-peptide complex

[unres.git] / ctest / matplotlib_fit_hist.py

diff --git a/ctest/matplotlib_fit_hist.py b/ctest/matplotlib_fit_hist.py
index 5a9bf37..754540c 100755 (executable)
--- a/ctest/matplotlib_fit_hist.py
+++ b/ctest/matplotlib_fit_hist.py
@@ -10,13 +10,16 @@ from math import sqrt
 import sys
 
 def prob_T(x,a):
-    g=np.float128(105.)
+    gg=np.float128(g)
+    aa=np.float128(10**(-gg-2)*a)
     Tr=np.float128(300.)
-    return a * (  x**((g-2)/2) * np.exp( -g*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[2])
 if (sys.argv[1] == '1L2Y_NH_GB000.stat' or sys.argv[1] == '1L2Y_NH_GB.stat'):
  x,y= np.loadtxt(sys.argv[1],usecols=(0,11),skiprows=10000,unpack=True)
 else: 
@@ -33,7 +36,7 @@ plt.xlabel('temperature')
 center = (bin[:-1] + bin[1:]) / 2
 #print bin
 #print center
-popt, pcov = curve_fit(prob_T, center, h, p0=10E-107)
+popt, pcov = curve_fit(prob_T, center, h)
 xfine = np.linspace(min(bin), max(bin), 100)
 #print popt