Contents
Fitting a 1-D PES using Python
Data set:
# Ar2 MP2 interaction energies
# Units: R : Angstrom
# E : cm^{-1}
#
# --------Using mid-bonds---
# R aDZ-noCP aTZ-noCP aQZ-noCP aDZ-CP aTZ-CP aQZ-CP
3.000 849.203 683.451 617.179 693.597 646.098 620.624
3.250 205.203 106.749 81.093 114.204 94.645 82.514
3.500 -14.944 -75.276 -82.840 -66.813 -76.199 -81.462
3.650 -63.415 -105.688 -109.144 -97.067 -104.068 -107.120
3.750 -78.828 -110.712 -112.906 -102.235 -108.261 -110.366
3.850 -85.975 -108.984 -110.509 -100.761 -106.027 -107.494
4.000 -87.113 -99.871 -100.789 -92.066 -96.351 -97.263
4.250 -76.570 -79.596 -79.416 -72.337 -75.132 -75.668
4.500 -61.796 -60.937 -59.607 -54.213 -55.757 -56.110
4.750 -48.339 -45.960 -44.074 -39.975 -40.728 -40.908
5.000 -37.453 -34.430 -32.644 -29.365 -29.747 -29.781
5.500 -22.077 -19.352 -18.364 -16.188 -16.262 -16.180
6.000 -12.663 -11.106 -10.662 -9.304 -9.278 -9.233import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt
unit_bohr2angstrom = 0.529177249
unit_hartree2cm = 219474.631371
def func_LJ(x,a,c):
# The Lennard-Jones model
return a/x**12 - c/x**6
def func_BM(x,a,b,c):
# The Born-Mayer model
#return a*np.exp(-b*x) - c/x**6
return np.exp(-b*(x-a)) - c/x**6
data = np.transpose(np.loadtxt('Ar2_MP2.txt'))
rvals = data[0]/unit_bohr2angstrom # convert from Angstrom to Bohr
evals = data[1]/unit_hartree2cm # convert from cm-1 to Hartree
# Standard deviations all set to 1.0
sigma = sigma = np.full(rvals.size,1.0)
print("sigma array: \n")
print(sigma)
# LJ fit
# =
# Starting values for parameters
x0=np.array([1.0,1.0])
# Do the fitting:
popt, pcov = opt.curve_fit(func_LJ,rvals,evals,p0=x0,sigma=sigma)
#
# popt : optimized parameters
# pcov : covariance matrix
print("\n\n")
print(" a = %12.4f c = %12.4f \n" % (popt[0],popt[1]))
print("Covariance matrix: \n")
ndim = pcov.shape[0]
for i in range(0,ndim):
for j in range(0,ndim):
pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
print(pcov)
print("\n\n")
#
# Plot
plt.plot(rvals,evals,'b-', label='data',marker='o')
plt.plot(rvals, func_LJ(rvals, *popt), 'g--', label='fit')
plt.show()
# BM fit
# =
# Starting values for parameters
x0=np.array([1.0,1.0,1.0])
# Do the fitting:
popt, pcov = opt.curve_fit(func_BM,rvals,evals,p0=x0,sigma=sigma)
#
# popt : optimized parameters
# pcov : covariance matrix
print("\n\n")
print(" a = %12.4f b = %12.4f c = %12.4f \n" % (popt[0],popt[1],popt[2]))
print("Covariance matrix: \n")
ndim = pcov.shape[0]
for i in range(0,ndim):
for j in range(0,ndim):
pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
print(pcov)
print("\n\n")
#
# Plot
plt.plot(rvals,evals,'b-', label='data',marker='o')
plt.plot(rvals, func_BM(rvals, *popt), 'g--', label='fit')
plt.show()
# Now try a weighting of exp(-(e-e0)/kT)
e0 = 80.0/unit_hartree2cm # both in units of evals
kT = 400.0/unit_hartree2cm
#sigma = np.exp(-np.abs(evals-e0)/kT)
sigma = np.exp(-np.abs(evals/e0))
print("sigma array: \n")
print(sigma)
# LJ fit
# =
# Starting values for parameters
x0=np.array([1.0,1.0])
# and re-fit
popt, pcov = opt.curve_fit(func_LJ,rvals,evals,p0=x0,sigma=sigma)
#
# popt : optimized parameters
# pcov : covariance matrix
print("\n\n")
print(" a = %12.4f c = %12.4f \n" % (popt[0],popt[1]))
print("Covariance matrix: \n")
ndim = pcov.shape[0]
for i in range(0,ndim):
for j in range(0,ndim):
pcov[i][j] = pcov[i][j]/np.sqrt(pcov[i][i]*pcov[j][j])
print(pcov)
print("\n\n")
#
# Plot
plt.plot(rvals,evals,'b-', label='data',marker='o')
plt.plot(rvals, func_LJ(rvals, *popt), 'g--', label='fit')
plt.show()
Covarianc
The fitting function **curve_fit** writes out not only the solution, but also the //covariance matrix//. This tells us how the parameters are correlated. To find out more about this have a look at this site. Here are a few examples of the //correlation matrix// which is defined in terms of the covariance matrix as follows: $$
- {\rm Corr}(x,y) = {\rm Cov}(x,y)/(\sigma_x \sigma_y)
$$ where the standard deviations are themselves part of the covariance matrix: $$
- \sigma_x = \sqrt{{\rm Cov}(x,x)}.
$$
Here's the correlation matrix for a fit to the Lennard-Jones function form: $$
V_{\rm LJ}(r) = \frac{A}{r{12}} - \frac{C}{r6}
$$
Correlation matrix: [[ 1. 0.9299438] [ 0.9299438 1. ]]
Ignore the $1.0$ terms on the diagonal. We are interested in the off-diagonal terms. For parameters that are not correlated these will be zero. Here we have off-diagonal terms nearly $+1$, so we conclude that the two parameters in the LJ function are heavily correlated, and will tend to change in the same way.
Conversely, here is the correlation matrix for the Born-Mayer function defined as: $$
- V_{\rm BM}(r) = \exp{(-\alpha(r-r_0))} - \frac{C}{r^6},
$$ with three free parameters: $\alpha$, $r_0$ and $C$.
Correlation matrix: [[ 1. 0.9904921 -0.87938874] [ 0.9904921 1. -0.93439446] [-0.87938874 -0.93439446 1. ]]
Here we see a mixed picture: $\alpha$ is heavily correlated with $r_0$. On the other hand, the correlation with $C$ is slightly weaker.