# Constructing the PES by interpolation#

The idea behind PES interpolation techniques is to construct the global PES of the molecular ground or excited-state from a set of data points (configurations) calculated quantum-mechanically (QM). For small molecules with few degrees of freedom, depending on the level of theory, the interpolated PES can have similar computational cost as constructing the PES ab initio. On the other hand, the use of an interpolated PES becomes very advantageous, especially when the molecule of interest is embedded in a larger environment (e.g. fluorescent probe in a protein), where this technique allows performing dynamics on the ground or excited state with quantum chemical accuracy, but at molecular mechanics computational cost [RP16].

The idea behind PES interpolation is to write the potential energy at a new configuration $$\mathbf{q}$$ as a weighted average between potentials defined from a set of ab initio data points $$\{\mathbf{q}_a\}$$:

$V(\mathbf{q}) = \sum_a w_a V_a(\mathbf{q})\,,$

where the weights $$w_a$$ depend on the distance ($$\boldsymbol{\Delta}_a$$) beween the new configuration $$\mathbf{q}$$ and a data configuration $$\mathbf{q}_a$$. The potentials $$V_a$$ are local harmonic potentials written in terms of Taylor expansions based on the energy ($$E_a$$), gradient ($$\mathbf{g}_a$$) and Hessian ($$\mathbf{H}_a$$):

$V_a(\mathbf{q}) = E_a + \boldsymbol{\Delta}^\mathrm{T}_a \mathbf{g}_a + \frac{1}{2} \boldsymbol{\Delta}^\mathrm{T}_a \mathbf{H}_a\boldsymbol{\Delta}_a\, .$

The vector $$\mathbf{q}$$ designating the position is in the internal coordinate system. The weights are provided by considering the distance between $$\mathbf{q}$$ and $$\mathbf{q}_a$$, such that the weight $$w_a$$ is larger for a data point $$\mathbf{q}_a$$ that is closer to $$\mathbf{q}$$. The weights $$w_a$$ are normalized weights:

$w_a = \frac{v_a}{\sum_i v_i}\,$

where the un-normalized weights $$v_a$$ can be calculated in different ways, for example a simple weighting scheme is:

$v_a = \frac{1}{|\mathbf{q}-\mathbf{q_a}|^{2p}}.$

The exponent $$2p$$ determines how fast the algorithm switches between data points [KR16].

Another approach to determine the un-normalized weights is the Shepard weighting scheme [KR16]:

$v_a = \left[ \left(\frac{|\mathbf{q}-\mathbf{q_a}|}{r_a}\right)^{2p} + \left(\frac{|\mathbf{q}-\mathbf{q_a}|}{r_a}\right)^{2q} \right]^{-1},$

where $$r_a$$ is the confidence radius for $$\mathbf{q_a}$$.

For a good overview of potential energy surfaces by interpolation, see [RP16].

In this tutorial we will get familiar with constructing the ground state potential energy surface (PES) of a molecule by interpolation.

## The ground state potential energy surface of water#

Let’s construct the ground state PES for the water molecule by interpolation. We need to carry out the following steps:

• Transform the coordinates, gradients and Hessians from Cartesian coordinates to internal coordinates,

• Calculate the ground state energy, gradient, and Hessian for a several molecular configurations which will serve as data points for interpolation,

• Determine the harmonic potentials and weights,

• Perform the interpolation and determine the potential energy at new configurations $$\mathbf{q}$$.

### Transforming from Cartesian to internal coordinates#

The steps required for the transformation of Cartesian coordinates to internal coordinates are described in details here. We will implement a set of routines to transform the gradient and Hessian to internal coordinates, making use of routines from geomeTRIC. Let’s start by setting up the molecule and calculating the energy, gradient and Hessian.

import veloxchem as vlx
import numpy as np
import sys
import py3Dmol as p3d
import geometric
from matplotlib import pyplot as plt

np.set_printoptions(suppress=True, precision=7)

* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.

* Warning * Setting MKL_THREADING_LAYER to "GNU".

molecule_string = """3
water
O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.895700000000       -0.316700000000
H    0.000000000000        0.000000000000        1.100000000000
"""
molecule = vlx.Molecule.from_xyz_string(molecule_string)
basis_set_label = 'sto-3g'

print("\nInitial Geometry of the water molecule:")
viewer = p3d.view(viewergrid=(1,1),width=300,height=200)
viewer.setStyle({'stick': {}})
viewer.rotate(-90,'y') # rotate the molecule to make it easier to see
viewer.show()

Initial Geometry of the water molecule:


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Calculating the desired properties:

scf_drv = vlx.ScfRestrictedDriver()

# Energy
scf_drv.compute(molecule, basis)
energy = scf_drv.get_scf_energy()

# Hessian
scf_hessian_drv = vlx.ScfHessianDriver(scf_drv)
scf_hessian_drv.compute(molecule, basis)
hessian = scf_hessian_drv.hessian


Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Restricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 8.6203186612 a.u.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.03 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -74.947250973879 a.u. Time: 0.04 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1       -74.947250981138    0.0000000000      0.00006046      0.00001095      0.00000000

                  2       -74.947250982569   -0.0000000014      0.00001394      0.00000251      0.00006130

                  3       -74.947250982655   -0.0000000001      0.00000159      0.00000036      0.00001805

                  4       -74.947250982655   -0.0000000000      0.00000011      0.00000002      0.00000116


*** SCF converged in 4 iterations. Time: 0.01 sec.


               Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy                       :      -74.9472509827 a.u.
Electronic Energy                  :      -83.5675696438 a.u.
Nuclear Repulsion Energy           :        8.6203186612 a.u.
------------------------------------

Ground State Information
------------------------
Charge of Molecule            :  0.0
Multiplicity (2S+1)           :  1.0
Magnetic Quantum Number (M_S) :  0.0

Spin Restricted Orbitals
------------------------

Molecular Orbital No.   1:
--------------------------
Occupation: 2.000 Energy:  -20.24331 a.u.
(   1 O   1s  :    -0.99)

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.23684 a.u.
(   1 O   1s  :    -0.24) (   1 O   2s  :     0.85) (   2 H   1s  :     0.18)

Molecular Orbital No.   3:
--------------------------
Occupation: 2.000 Energy:   -0.59041 a.u.
(   1 O   1p-1:    -0.41) (   1 O   1p0 :     0.44) (   2 H   1s  :    -0.45)
(   3 H   1s  :     0.43)

Molecular Orbital No.   4:
--------------------------
Occupation: 2.000 Energy:   -0.43151 a.u.
(   1 O   2s  :    -0.49) (   1 O   1p-1:     0.60) (   1 O   1p0 :     0.49)
(   2 H   1s  :     0.24) (   3 H   1s  :     0.35)

Molecular Orbital No.   5:
--------------------------
Occupation: 2.000 Energy:   -0.38490 a.u.
(   1 O   1p+1:     1.00)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.000 Energy:    0.50463 a.u.
(   1 O   2s  :    -0.69) (   1 O   1p-1:    -0.34) (   1 O   1p0 :    -0.66)
(   2 H   1s  :     0.44) (   3 H   1s  :     0.94)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.000 Energy:    0.71372 a.u.
(   1 O   2s  :    -0.41) (   1 O   1p-1:    -0.76) (   1 O   1p0 :     0.58)
(   2 H   1s  :     1.08) (   3 H   1s  :    -0.41)



=====================


                                              Molecular Geometry (Angstroms)
================================

Atom         Coordinate X          Coordinate Y          Coordinate Z

O           0.000000000000        0.000000000000        0.000000000000
H           0.000000000000        0.895700000000       -0.316700000000
H           0.000000000000        0.000000000000        1.100000000000

-------------------------

O           0.000000000000        0.085297089619       -0.100972944693
H           0.000000000014       -0.061749632536       -0.007085216623
H          -0.000000000014       -0.023547457801        0.108058171904


                                   *** Time spent in gradient calculation: 1.29 sec ***



SCF Hessian Driver
====================


                                   *** Time spent in Hessian calculation: 13.68 sec ***



Now that we have an energy, gradient and Hessian, let’s define the Z-matrix for the water molecule and write the routines required to transform from Cartesian coordinates to the internal coordinates defined by the Z-matrix. We will use geomeTRIC to set up the internal coordinates.

#            O-H1    O-H2   H1-O-H2
z_matrix = [[0, 1], [0,2], [1,0,2]]

internal_coordinates = []

for z in z_matrix:
if len(z) == 2:
# define a bond distance object
q = geometric.internal.Distance(z[0], z[1])

elif len(z) == 3:
# define a bond angle object
q = geometric.internal.Angle(z[0], z[1], z[2])

else:
# define a dihedral angle object
q = geometric.internal.Dihedral(z[0], z[1], z[2], z[3])

internal_coordinates.append(q)


The internal coordinates defined through geomeTRIC have routines that allow us to calculate derivatives with respect to Cartesian coordinates ($$\mathbf{x}$$), as required for constructing the Wilson $$\mathbf{B}$$-matrix:

$B_{ij} = \frac{\partial q_i}{\partial x_j}\,.$

Let us use these routines to get the $$\mathbf{B}$$-matrix.

def get_b_matrix(molecule, internal_coordinates, z_matrix):
n_atoms = molecule.number_of_atoms()

# number of Cartesian coordinates
n_cart = 3 * n_atoms

# number of internal coordinates
n_int = len(internal_coordinates)

# initialize the B-matrix
b_matrix = np.zeros((n_int, n_cart))

# Cartesian coordinates
coords = molecule.get_coordinates().reshape( 3 * n_atoms)

# calculate the derivatives of q with respect to coords
i = 0 # i runs over all the internal coordinates
for q in internal_coordinates:
deriv_q = q.derivative(coords)
# add the derivative values in the right spots of the b_matrix;
# we need the atom indices from the Z-matrix.
for a in z_matrix[i]:
# from 3 * atom_index to 3 * atom_index + 3
b_matrix[i, 3 * a : 3 * a + 3 ] = deriv_q[a]
i += 1
return b_matrix

b_mat = get_b_matrix(molecule, internal_coordinates, z_matrix)


In order to transform from Cartesian gradient to internal gradient, we would need the inverse of $$\mathbf{B}$$:

$\frac{\partial V}{\partial q_i} = \sum_j \frac{\partial V}{\partial x_j} \frac{\partial x_j}{\partial q_i} = \sum_j \frac{\partial V}{\partial x_j} B_{ji}^{-},$

with $$\mathbf{BB}^{-}=\mathbf{I}$$. Because $$\mathbf{B}$$ is a rectangular matrix, we need to obtain its inverse carefully. From

$\mathbf{B}\mathbf{B}^\mathrm{T} \left( \mathbf{B}\mathbf{B}^\mathrm{T} \right)^{-} = \mathbf{I},$

we can easily see that $$\mathbf{B}^{-} = \mathbf{B}^\mathrm{T} \left( \mathbf{B}\mathbf{B}^\mathrm{T} \right)^{-}$$ holds, so to get the inverse we have to first construct a square matrix $$\mathbf{G}=\mathbf{B}\mathbf{B}^\mathrm{T}$$:

g_mat = ...

g_mat = np.matmul(b_mat, b_mat.T)


$$\mathbf{G}$$ is diagnoalized to get its eigenvalues and eigenvectors. The inverse of the non-zero eigenvalues and corresponding eigenvectors will be then used to construct the $$\mathbf{G}^{-}$$ matrix.

g_val, g_vec = ...

# Invert the non-zero eigenvalues
...
...

g_minus = np.linalg.multi_dot([g_vec, g_val_inverse_mat, g_vec.T])

g_val, g_vec = np.linalg.eig(g_mat)
g_val_inverse = []
for g in g_val:
if abs(g) > 1e-7:
g_val_inverse.append(1.0 / g)
else:
g_val_inverse.append(0.0)
g_val_inverse_mat = np.diag(np.array(g_val_inverse))
g_minus = np.linalg.multi_dot([g_vec, g_val_inverse_mat, g_vec.T])


Let’s make this piece of code a general routine that we can use for an arbitrary $$\mathbf{B}$$ matrix.

def get_g_minus(b_matrix):

...

return g_minus

def get_g_minus(b_matrix):
g_mat = np.matmul(b_matrix, b_matrix.T)
g_val, g_vec = np.linalg.eig(g_mat)
g_val_inverse = []
for g in g_val:
if abs(g) > 1e-7:
g_val_inverse.append(1.0 / g)
else:
g_val_inverse.append(0.0)

g_val_inverse_mat = np.diag(np.array(g_val_inverse))
g_minus = np.linalg.multi_dot([g_vec, g_val_inverse_mat, g_vec.T])

return g_minus

g = get_g_minus(b_mat)


Now we have all the ingredients to transform the gradient, but we need the second order derivatives of $$\mathbf{q}$$ with respect to $$\mathbf{x}$$ to be able to also transform the Hessian. Let’s write a routine to calculate the second-order B-matrix, $$\mathbf{B}^{(2)}$$.

def get_b2_matrix(molecule, internal_coordinates, z_matrix):

# Use as an example the routine to construct the B-matrix
# and write a similar routine for the matrix of second order
# derivatives

return b2_matrix

def get_b2_matrix(molecule, internal_coordinates, z_matrix):
n_atoms = molecule.number_of_atoms()

# number of Cartesian coordinates
n_cart = 3 * n_atoms

# number of internal coordinates
n_int = len(internal_coordinates)

# initialize the B-matrix
b2_matrix = np.zeros((n_int, n_cart, n_cart))

# Cartesian coordinates
coords = molecule.get_coordinates().reshape( 3 * n_atoms)

# calculate the derivatives of q with respect to coords
i = 0 # i runs over all the internal coordinates
for q in internal_coordinates:
second_deriv_q = q.second_derivative(coords)
# add the derivative values in the right spots of the b_matrix;
# we need the atom indices from the Z-matrix.
for a in z_matrix[i]:
for b in z_matrix[i]:
# from 3 * atom_index_a to 3 * atom_index_a + 3
# and from 3 * atom_index_b to 3 * atom_index_b + 3
b2_matrix[i, 3 * a : 3 * a + 3,
3 * b : 3 * b + 3] = second_deriv_q[a, :, b, :]
i += 1
return b2_matrix

b2_matrix = get_b2_matrix(molecule, internal_coordinates, z_matrix)


Let’s write the routines which transform the gradient and the Hessian from Cartesian coordinates to internal coordinates. We have to implement the following equations:

$\mathbf{g}_q = \mathbf{G}^{-}\mathbf{B}\mathbf{g}_x$
$\mathbf{H}_q = \mathbf{G}^{-}\mathbf{B}\big[ \mathbf{H}_x - \mathbf{g}_q\mathbf{B}^{(2)}\big]\mathbf{B}^\mathrm{T}\mathbf{G}^{-\mathrm{T}}$
def convert_gradient_to_internal(gradient, b_matrix, g_minus):

# Implement the equation to convert the gradient from Cartesian
# to internal coordinates.


def convert_gradient_to_internal(gradient, b_matrix, g_minus):

# reshape the Cartesian gradient to match the dimensions of the
# g_minus and b_matrix

# convert to internal coordinates

# To check if the gradient in internal coordinates is correct, you can calculate selected elements numerically.


Note

Has the conversion worked correctly? How can you check?

Now, let’s write a routine to convert the Hessian.

def convert_Hessian_to_internal(hessian, b_matrix, g_minus, b2_matrix, grad_q):
# calculate term related to the second order B-matrix

# convert the Hessian from Cartesian to internal coordinates
hessian_q = ...

return hessian_q

def convert_Hessian_to_internal(hessian, b_matrix, g_minus, b2_matrix, grad_q):
# calculate term related to the second order B-matrix

# convert the Hessian from Cartesian to internal coordinates
hessian_q = np.linalg.multi_dot([g_minus, b_matrix, hessian - gradq_b2mat, b_matrix.T, g_minus.T])
return hessian_q

hessian_q = convert_Hessian_to_internal(hessian, b_mat, g_minus, b2_matrix, grad_q)

# To check if the conversion of the Hessian works correctly
# you can calculate selected Hessian matrix elements numerically.


### Ground-state PES by interpolation#

Now that we have set up the routines for conversion between Cartesian and internal coordiantes, we can start collecting data points for interpolation. Let’s start with stretching one of the OH bonds and calculate three new points: one when the bond is contracted, one around equilibrium bond length, and one when the bond is stretched.

mol_template = """
O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.895700000000       -0.316700000000
H    0.000000000000        0.000000000000        OHdist
"""

# list of bond distances
distlist = [0.75, 0.95, 1.35]

# empty list to save the energies, gradients and hessians
energies = []
cart_hessians = []

# list for molecular configurations
molecules = []

for oh in distlist:
print("Calculating the energies, gradients and Hessians for OH = ", oh, "A...")

# Create new molecule
mol_str = mol_template.replace("OHdist", str(oh))

# set-up an scf driver
new_scf_drv = vlx.ScfRestrictedDriver()
new_scf_drv.ostream.state = False # To disable printout
new_scf_drv.compute(new_molecule, basis)

# calculate the Hessian
new_scf_hessian_drv = vlx.ScfHessianDriver(new_scf_drv)
new_scf_hessian_drv.ostream.state = False
new_scf_hessian_drv.compute(new_molecule, basis)

# save the results:
energy = new_scf_drv.get_scf_energy()
cart_hessian = new_scf_hessian_drv.hessian

energies.append(energy)
cart_hessians.append(cart_hessian)
molecules.append(new_molecule)

Calculating the energies, gradients and Hessians for OH =  0.75 A...

Calculating the energies, gradients and Hessians for OH =  0.95 A...

Calculating the energies, gradients and Hessians for OH =  1.35 A...


Now we have to transform the gradient and Hessian from Cartesian to internal coordinates in order to be able to use them for interpolation.

# number of data points:
n_data_pts = len(energies)

# empty lists for the internal gradients and Hessians:
int_hessians = []

for i in range(n_data_pts):
...
# convert Hessians
...

# number of data points:
n_data_pts = len(energies)

# empty lists for the internal gradients and Hessians:
int_hessians = []

for i in range(n_data_pts):
hess_x = cart_hessians[i]
mol = molecules[i]
b_mat = get_b_matrix(mol, internal_coordinates, z_matrix)
g_minus = get_g_minus(b_mat)
b2_mat = get_b2_matrix(mol, internal_coordinates, z_matrix)

hess_q = convert_Hessian_to_internal(hess_x, b_mat, g_minus, b2_mat, grad_q)
int_hessians.append(hess_q)


Now that we have the gradients and Hessians in internal coordinates, let’s write a routine which constructs the harmonic potential, given a new molecular configuration $$\mathbf{q}$$.

$V_a(\mathbf{q}) = E_a + \boldsymbol{\Delta}^\mathrm{T}_a \mathbf{g}_a + \frac{1}{2} \boldsymbol{\Delta}^\mathrm{T}_a \mathbf{H}_a\boldsymbol{\Delta}_a\, .$
def get_harmonic_potential(q_a, e_a, g_a, h_a, new_q):
# q_a are the internal coordinates where the energy, gradient and
# hessian have been calculated;
dist = ...
pot = ...
return pot

def get_harmonic_potential(q_a, e_a, g_a, h_a, new_q):
# coords_a are the coordinates where the energy, gradient and
# hessian have been calculated;
dist = new_q - q_a
pot = e_a + np.matmul(dist.T, g_a) + 0.5 * np.linalg.multi_dot([dist.T, h_a, dist])
return pot


To be able to use the routine, we first have to create the array of internal coordinates $$\mathbf{q}_a$$. Let’s add a routine to calculate the values of the internal coordinates we defined previously using geomeTRIC objects. We will need the Cartesian coordinates of the molecule in the different configurations.

Note

The Cartesian and internal coordinates used here are in atomic units!

def get_internal_coordinates_array(internal_coordinates, cart_coords):
q_list = []
for q in internal_coordinates:
q_list.append(q.value(cart_coords))
return np.array(q_list)

# list of internal coordinates with our interpolation data points:
qs = []

for molec in molecules:
cart_coords = molec.get_coordinates()
q_a = get_internal_coordinates_array(internal_coordinates, cart_coords)
qs.append(q_a)

def get_internal_coordinates_array(internal_coordinates, cart_coords):
q_list = []
for q in internal_coordinates:
q_list.append(q.value(cart_coords))
return np.array(q_list)

# list of internal coordinates with our interpolation data points:
qs = []

for molec in molecules:
cart_coords = molec.get_coordinates()
q_a = get_internal_coordinates_array(internal_coordinates, cart_coords)
qs.append(q_a)


Let’s construct the harmonic potentials based on the three data points we calculated previously. To plot the results let’s use an array of new O-H bond distances.

distlist = np.arange(0.6,2.0,0.05)

# list of harmonic potentials
harmonic_potentials = []

for i in range(n_data_pts):
q_a = qs[i]
e_a = energies[i]
h_a = int_hessians[i]

potential = []

for oh in distlist:
#print("Calculating the harmonic potential for... ", oh)

# Create new molecule
mol_str = mol_template.replace("OHdist", str(oh))
cart_coords = new_molecule.get_coordinates()
q = get_internal_coordinates_array(internal_coordinates, cart_coords)

# Calculate harmonic potential at point q
pot = get_harmonic_potential(q_a, e_a, g_a, h_a, q)
potential.append(pot)
harmonic_potentials.append(np.array(potential))


Let’s plot the results.

plt.figure(figsize=(6,4))
plt.title('Local Harmonic Potentials for Interpolation')
plt.plot(distlist, harmonic_potentials[0], label="0.75 Å")
plt.plot(distlist, harmonic_potentials[1], label="0.95 Å")
plt.plot(distlist, harmonic_potentials[2], label="1.35 Å")

plt.xlabel('O-H bond length (Å)')
plt.ylabel('Energy (H)')

plt.legend()
plt.show()


Let’s use these for interpolation and compare to the PES calculated with HF. We need to implement the following equations:

$V(\mathbf{q}) = \sum_a w_a V_a(\mathbf{q})\,,$
$w_a = \frac{v_a}{\sum_i v_i}\,$
$v_a = \frac{1}{|\mathbf{q}-\mathbf{q_a}|^{2p}}.$
def get_interpolated_pes(new_point, data_points, energies, gradients, hessians, exponent):
# sum of all interpolation weights -- to normalize them
sum_weights = 0
weights = []
potentials = []

n_points = len(data_points)

for i in range(n_points):
# Calculate the norm of the distance vector

# Calculate the un-normalized weight

# Calculate the harmonic potential

# Calculate normalized weights

# Calculate interpolated energy at new_point

return new_energy

def get_interpolated_pes(new_point, data_points, energies, gradients, hessians, exponent):
# sum of all interpolation weights -- to normalize them
sum_weights = 0
weights = []
potentials = []

n_points = len(data_points)

for i in range(n_points):
q_a = data_points[i]
e_a = energies[i]
h_a = hessians[i]
dist = np.linalg.norm(new_point - q_a)
weight = 1 / dist**exponent
pot = get_harmonic_potential(q_a, e_a, g_a, h_a, new_point)
#print(pot)
sum_weights += weight
weights.append(weight)
potentials.append(pot)

new_energy = 0
for i in range(len(data_points)):
new_energy += weights[i] / sum_weights * potentials[i]
#print(en)

return new_energy


Let’s see how well the interpolated energies compare to HF energies.

scf_energies = []
im_energies = []

for oh in distlist:

# Create new molecule
...

# set-up an scf driver and calculate the energy
...

# calculate energy by interpolation
...

scf_energies = []
im_energies = []

for oh in distlist:

# Create new molecule
mol_str = mol_template.replace("OHdist", str(oh))

# set-up an xtb driver, gradient driver and hessian driver
new_scf_drv = vlx.ScfRestrictedDriver()
new_scf_drv.ostream.state = False
new_scf_drv.compute(new_molecule, basis)

# save the results:
energy = new_scf_drv.get_scf_energy()

scf_energies.append(energy)

# calculate energies by interpolation
coords = new_molecule.get_coordinates()
new_q = get_internal_coordinates_array(internal_coordinates, coords)
im_energy = get_interpolated_pes(new_q, qs, energies, int_gradients, int_hessians, 12)

im_energies.append(im_energy)

qm_points = [0.75, 0.95, 1.35]

plt.figure(figsize=(6,4))
plt.title('Energy during O-H dissociation')
plt.plot(distlist, scf_energies, label='HF')
plt.plot(distlist, im_energies, label='Interpolated')
plt.plot(qm_points, energies, 'o', label='QM data points',  color="dimgrey")

plt.xlabel('O-H bond length (Å)')
plt.ylabel('Energy (H)')

plt.legend()
plt.tight_layout(); plt.show()

plt.figure(figsize=(6,4))
plt.title('Interpolated Energy during O-H dissociation')
plt.plot(distlist, im_energies, label='IM')
plt.plot(distlist, harmonic_potentials[0], label="0.75 Å")
plt.plot(distlist, harmonic_potentials[1], label="0.95 Å")
plt.plot(distlist, harmonic_potentials[2], label="1.35 Å")
plt.plot(qm_points, energies, 'o', label='QM data points',  color="dimgrey")
plt.axis(ymin=-75.1, ymax=-74.5)

plt.xlabel('O-H bond length (Å)')
plt.ylabel('Energy (H)')

plt.legend()
plt.tight_layout(); plt.show()


Note

• What happens if you calculate the energy of a configuration which is outside the range of QM data points (eg. an OH bond length larger than the largest value you included)?

• What happens if you include more or fewer data points in your data base?