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     Tr=np.float128(300.)
  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) ) )
  18
  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)
  22 g=int(sys.argv[1])
  23 with open('file_GB000.stat','r') as f:
  24   for line in f:
  25     pass
  26   ncolumns=len(line.split())
  27
  28 if ncolumns==14:
  29  x,y= np.loadtxt('file_GB000.stat',usecols=(0,10),skiprows=10,unpack=True)
  30  x1,e,r= np.loadtxt('file_GB000.stat',usecols=(0,3,5),skiprows=0,unpack=True)
  31 else:
  32  x,y= np.loadtxt('file_GB000.stat',usecols=(0,6),skiprows=10,unpack=True)
  33  x1,e= np.loadtxt('file_GB000.stat',usecols=(0,3),skiprows=0,unpack=True)
  34
  35 h,bin=np.histogram(y,bins=50,density=True)
  36
  37 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
  38 plt.xlim(min(bin), max(bin))
  39 plt.ylim(0,max(h))
  40 plt.ylabel('probality')
  41 plt.xlabel('temperature')
  42
  43 center = (bin[:-1] + bin[1:]) / 2
  44 #print bin
  45 #print center
  46 popt, pcov = curve_fit(prob_T, center, h)
  47 xfine = np.linspace(min(bin), max(bin), 100)
  48 #print popt
  49
  50 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
  51 print '%15.10f' % chi_squared
  52
  53 #print xfine
  54 #print prob_T(xfine,popt[0])
  55 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
  56 plt.savefig('temp_hist.png')
  57 #plt.show()
  58
  59 plt.clf()
  60 plt.xlabel('step')
  61 plt.ylabel('potential energy')
  62 plt.plot(x1,e,'.')
  63 plt.savefig('md_ene.png')
  64
  65 if ncolumns==14:
  66  plt.clf()
  67  plt.xlabel('step')
  68  plt.ylabel('RMSD')
  69  plt.plot(x1,r,'.')
  70  plt.savefig('md_rms.png')
  71