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#

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)
Hide code cell output
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 4.
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.read_xyz_string(molecule_string)
basis_set_label = 'sto-3g'
basis = vlx.MolecularBasis.read(molecule, basis_set_label)
* Info * Reading basis set from file: /home/thomas/miniconda3/envs/echem/lib/python3.10/site-packages/veloxchem/basis/STO-3G
                                                                                                                          
                                              Molecular Basis (Atomic Basis)                                              
                                             ================================                                             
                                                                                                                          
                               Basis: STO-3G                                                                              
                                                                                                                          
                               Atom Contracted GTOs           Primitive GTOs                                              
                                                                                                                          
                                O   (2S,1P)                   (6S,3P)                                                     
                                H   (1S)                      (3S)                                                        
                                                                                                                          
                               Contracted Basis Functions : 7                                                             
                               Primitive Basis Functions  : 21                                                            
                                                                                                                          
Hide code cell source
print("\nInitial Geometry of the water molecule:")
viewer = p3d.view(viewergrid=(1,1),width=300,height=200)
viewer.addModel(molecule_string, 'xyz', viewer=(0,0))
viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({'stick': {}, 'sphere': {'scale':0.25}})
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_results= scf_drv.compute(molecule, basis)
energy = scf_drv.get_scf_energy()

# Gradient
scf_grad_drv = vlx.ScfGradientDriver()
scf_grad_drv.compute(molecule, basis, scf_drv)
gradient = scf_grad_drv.gradient

# Hessian
scf_hessian_drv = vlx.ScfHessianDriver()
scf_hessian_drv.compute(molecule, basis, scf_drv)
hessian = scf_hessian_drv.hessian
Hide code cell output
                                                                                                                          
                                            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.01 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.947250981137    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.02 sec.                                                         
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :      -74.9472509827 a.u.                                              
               Electronic Energy                  :      -83.5675696438 a.u.                                              
               Nuclear Repulsion Energy           :        8.6203186612 a.u.                                              
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000001106 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)                                                        
                                                                                                                          
                                                                                                                          
                                                   SCF Gradient Driver                                                    
                                                  =====================                                                   
                                                                                                                          
                               Gradient Type               : Numerical                                                    
                               Numerical Method            : Symmetric Difference Quotient                                
                               Finite Difference Step Size : 0.001 a.u.                                                   
                                                                                                                          
                                              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                         
                                                                                                                          
                                            Numerical Gradient (Hartree/Bohr)                                             
                                           -----------------------------------                                            
                                                                                                                          
                          Atom           Gradient X            Gradient Y            Gradient Z                           
                                                                                                                          
                           O           0.000000000000        0.085297089612       -0.100972944672                         
                           H           0.000000000043       -0.061749632522       -0.007085216652                         
                           H           0.000000000000       -0.023547457772        0.108058171861                         
                                                                                                                          
                                   *** Time spent in gradient calculation: 1.69 sec ***                                   
                                                                                                                          
                                                                                                                          
                                                    SCF Hessian Driver                                                    
                                                   ====================                                                   
                                                                                                                          
                                   *** Time spent in Hessian calculation: 15.81 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 = ...
Hide code cell source
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])
Hide code cell source
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
Hide code cell source
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
Hide code cell source
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.
    
    return gradient_q

Hide code cell source
def convert_gradient_to_internal(gradient, b_matrix, g_minus):
    n_atoms = gradient.shape[0]
    
    # reshape the Cartesian gradient to match the dimensions of the 
    # g_minus and b_matrix
    grad_x = gradient.reshape((3 * n_atoms))
    
    # convert to internal coordinates
    grad_q = np.linalg.multi_dot([g_minus, b_matrix, grad_x])
    return grad_q

grad_q = convert_gradient_to_internal(gradient, b_mat, g_minus)

# 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
    gradq_b2mat = ...
    
    # convert the Hessian from Cartesian to internal coordinates
    hessian_q = ...
    
    return hessian_q
Hide code cell source
def convert_Hessian_to_internal(hessian, b_matrix, g_minus, b2_matrix, grad_q):
    # calculate term related to the second order B-matrix
    gradq_b2mat = np.einsum("i,ixy->xy", grad_q, b2_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_gradients = []
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))
    new_molecule = vlx.Molecule.read_molecule_string(mol_str, units='angstrom')
        
    # set-up an scf driver
    new_scf_drv = vlx.ScfRestrictedDriver()
    new_scf_drv.ostream.state = False # To disable printout
    new_scf_results = new_scf_drv.compute(new_molecule, basis)
    
    # calculate the gradient
    new_scf_grad_drv = vlx.ScfGradientDriver()
    new_scf_grad_drv.ostream.state = False
    new_scf_grad_drv.compute(new_molecule, basis, new_scf_drv)
    
    # calculate the Hessian
    new_scf_hessian_drv = vlx.ScfHessianDriver()
    new_scf_hessian_drv.ostream.state = False
    new_scf_hessian_drv.compute(new_molecule, basis, new_scf_drv)
    
    # save the results:
    energy = new_scf_drv.get_scf_energy()
    cart_gradient = new_scf_grad_drv.gradient
    cart_hessian = new_scf_hessian_drv.hessian
    
    energies.append(energy)
    cart_gradients.append(cart_gradient)
    cart_hessians.append(cart_hessian)
    molecules.append(new_molecule) 
Hide code cell output
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_gradients = []
int_hessians = []

for i in range(n_data_pts):
    # convert gradients
    ...
    # convert Hessians
    ...
Hide code cell source
# number of data points:
n_data_pts = len(energies)

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

for i in range(n_data_pts):
    grad_x = cart_gradients[i]
    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)
    
    grad_q = convert_gradient_to_internal(grad_x, b_mat, g_minus)
    hess_q = convert_Hessian_to_internal(hess_x, b_mat, g_minus, b2_mat, grad_q)
    int_gradients.append(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
Hide code cell source
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)
Hide code cell source
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]
    g_a = int_gradients[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))
        new_molecule = vlx.Molecule.read_molecule_string(mol_str, units='angstrom')
        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()
Hide code cell output
../../_images/9c4837a407eb3fb47dc1ade17cc5a2c4ac083104d0e34701e33d250125bd481f.png

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
Hide code cell source
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]
        g_a = gradients[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
    ...
Hide code cell source
scf_energies = []
im_energies = []

for oh in distlist:
    
    # Create new molecule
    mol_str = mol_template.replace("OHdist", str(oh))
    new_molecule = vlx.Molecule.read_molecule_string(mol_str, units='angstrom')
    
    # set-up an xtb driver, gradient driver and hessian driver
    new_scf_drv = vlx.ScfRestrictedDriver()
    new_scf_drv.ostream.state = False
    new_scf_results = 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()
Hide code cell output
../../_images/4db05dab3a710b4c80b21751eacea1660473bb452c9940d149972a218c2ea7b5.png
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()
Hide code cell output
../../_images/4d16924d01fcfade87066b9cae69777749c6d405c4297c83898075e0bf5f6f64.png

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?