attachment:fit_curve_v1.py of AJMPublic/teaching/electronic-structure/fitting-1d-pes-python

   1 #!/usr/bin/python3.5
   2 
   3 import math
   4 import matplotlib.pyplot as plt
   5 import numpy as np
   6 import scipy.optimize as opt
   7 
   8 unit_bohr2angstrom = 0.529177249
   9 unit_hartree2cm    = 219474.631371
  10 
  11 def func_LJ(x,a,c):
  12     # The Lennard-Jones model
  13     return a/x**12 - c/x**6
  14 
  15 def func_BM(x,a,b,c):
  16     # The Born-Mayer model
  17     #return a*np.exp(-b*x) - c/x**6
  18     return np.exp(-b*(x-a)) - c/x**6
  19 
  20 data = np.transpose(np.loadtxt('Ar2_MP2.txt'))
  21 rvals = data[0]/unit_bohr2angstrom # convert from Angstrom to Bohr
  22 evals = data[1]/unit_hartree2cm    # convert from cm-1 to Hartree
  23 
  24 
  25 # Standard deviations all set to 1.0
  26 sigma = sigma = np.full(rvals.size,1.0)
  27 
  28 print("sigma array: \n")
  29 print(sigma)
  30 
  31 # LJ fit
  32 # ======
  33 # Starting values for parameters
  34 x0=np.array([1.0,1.0])
  35 # Do the fitting:
  36 popt, pcov = opt.curve_fit(func_LJ,rvals,evals,p0=x0,sigma=sigma)
  37 #
  38 # popt : optimized parameters
  39 # pcov : covariance matrix
  40 print("\n\n")
  41 print(" a = %12.4f  c = %12.4f  \n" % (popt[0],popt[1]))
  42 print("Covariance matrix: \n")
  43 ndim = pcov.shape[0]
  44 for i in range(0,ndim):
  45     for j in range(0,ndim):
  46         pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
  47 print(pcov)
  48 print("\n\n")
  49 #
  50 # Plot
  51 plt.plot(rvals,evals,'b-', label='data',marker='o')
  52 plt.plot(rvals, func_LJ(rvals, *popt), 'g--', label='fit')
  53 plt.show()
  54 
  55 # BM fit
  56 # ======
  57 # Starting values for parameters
  58 x0=np.array([1.0,1.0,1.0])
  59 # Do the fitting:
  60 popt, pcov = opt.curve_fit(func_BM,rvals,evals,p0=x0,sigma=sigma)
  61 #
  62 # popt : optimized parameters
  63 # pcov : covariance matrix
  64 print("\n\n")
  65 print(" a = %12.4f  b = %12.4f  c = %12.4f  \n" % (popt[0],popt[1],popt[2]))
  66 print("Covariance matrix: \n")
  67 ndim = pcov.shape[0]
  68 for i in range(0,ndim):
  69     for j in range(0,ndim):
  70         pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
  71 print(pcov)
  72 print("\n\n")
  73 #
  74 # Plot
  75 plt.plot(rvals,evals,'b-', label='data',marker='o')
  76 plt.plot(rvals, func_BM(rvals, *popt), 'g--', label='fit')
  77 plt.show()
  78 
  79 
  80 # Now try a weighting of exp(-(e-e0)/kT)
  81 e0 =   80.0/unit_hartree2cm  # both in units of evals
  82 kT =  400.0/unit_hartree2cm
  83 #sigma = np.exp(-np.abs(evals-e0)/kT)
  84 sigma = np.exp(-np.abs(evals/e0))
  85 
  86 print("sigma array: \n")
  87 print(sigma)
  88 
  89 # LJ fit
  90 # ======
  91 # Starting values for parameters
  92 x0=np.array([1.0,1.0])
  93 # and re-fit
  94 popt, pcov = opt.curve_fit(func_LJ,rvals,evals,p0=x0,sigma=sigma)
  95 #
  96 # popt : optimized parameters
  97 # pcov : covariance matrix
  98 print("\n\n")
  99 print(" a = %12.4f  c = %12.4f  \n" % (popt[0],popt[1]))
 100 print("Covariance matrix: \n")
 101 ndim = pcov.shape[0]
 102 for i in range(0,ndim):
 103     for j in range(0,ndim):
 104         pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
 105 print(pcov)
 106 print("\n\n")
 107 #
 108 # Plot
 109 plt.plot(rvals,evals,'b-', label='data',marker='o')
 110 plt.plot(rvals, func_LJ(rvals, *popt), 'g--', label='fit')
 111 plt.show()