Overlap Integrals

We can now proceed to calculate molecular integrals, starting with overlap integrals. An overlap integral between two centers $A$ and $B$ are defined as $$ S_{ijk,lmn}^{A,B} = \int g_{ijk}(\vec{r};\alpha, \vec{A}) \ g_{lmn}(\vec{r};\beta, \vec{B}) \ \mathrm{d}^3 r $$ Because Cartesian Gaussians are factorizable, we can calculate the overlap for each Cartesian direction, i.e., $$ S_{ijk,lmn}^{A,B} = S_{il}^{A_x, B_x}\ S_{jm}^{A_y, B_y}\ S_{kn}^{A_z, B_z} $$ Therefore, we only have to deal with overlaps between 1D Gaussians.

Hermite Gaussian Overlaps

We now examine the overlap between two (unnormalized) 1D Hermite Gaussian: $$ \begin{align} E_{k l}^{A_x, B_x} &= \int h_{k}(x;\alpha, A_x) \ h_{l}(x;\beta, B_x) \ \mathrm{d} x \\ &= \int \left( \frac{\partial}{\partial A_x} \right)^k \mathrm{e}^{-\alpha(x-A_x)^2} \left( \frac{\partial}{\partial B_x} \right)^l \mathrm{e}^{-\beta(x-B_x)^2} \mathrm{d} x \\ &= \left( \frac{\partial}{\partial A_x} \right)^k \left( \frac{\partial}{\partial B_x} \right)^l \int \mathrm{e}^{-\alpha(x-A_x)^2} \mathrm{e}^{-\beta(x-B_x)^2} \mathrm{d} x \\ &= \left( \frac{\partial}{\partial A_x} \right)^k \left( \frac{\partial}{\partial B_x} \right)^l E_{00}^{A_x,B_x} \end{align} $$ By using the Gaussian product theorem the integral $E_{00}^{A_x,B_x}$ is reduced to $$ \begin{align} E_{00}^{A_x,B_x} &= \mathrm{e}^{-\mu X_{AB}^2} \int \mathrm{e}^{-p (x - P_x)^2} \\ &= \mathrm{e}^{-\mu X_{AB}^2} \sqrt{\frac{\pi}{p}} \end{align} $$ with $$ \begin{align} p &= \alpha + \beta \\ \mu &= \frac{\alpha \beta}{\alpha + \beta} \\ X_{AB} &= A_x - B_x \\ P &= \frac{\alpha A_x + \beta B_x}{\alpha + \beta} \end{align} $$ Substituting the variables back, we get $$ E_{00}^{A_x,B_x} = \sqrt{\frac{\pi}{\alpha + \beta}}\ \exp \left( -\frac{\alpha \beta}{\alpha + \beta} (A_x - B_x)^2 \right) $$ Differentiate this expression with respect to $A_x$ and $B_x$ will deliver us with all possible overlap integrals between Hermite Gaussians.

Cartesian Gaussian Overlaps

By expanding Cartesian Gaussians into Hermite Gaussians, we can easily obtain their overlaps, i.e., $$ \begin{align} S_{i j}^{A_x, B_x} &= \int g_i(x; \alpha, A_x)\ g_j(x; \beta, B_x)\ \mathrm{d}x \\ &= \int \sum_{k} c_{ki} h_k(x; \alpha, A_x) \sum_{l} c_{lj} h_l(x; \beta, B_x) \mathrm{d}x \\ &= \sum_{k} \sum_{l} c_{ki} c_{lj} \int h_k(x; \alpha, A_x)\ h_l(x; \beta, B_x)\ \mathrm{d}x \\ &= \sum_{k} \sum_{l} c_{ki} c_{lj} E_{kl}^{A_x,B_x} \end{align} $$ So, we can obtain the Cartesian Gaussian overlaps by linearly combining Hermite Gaussian Overlaps with the corresponding Hermite expansion coefficients.

We shall now use SymPy to generate formulas for cartesian overlaps.

Code Generation

We start by importing necessary modules, including our function for calculating Hermite expansion coefficients. It is assumed here that this function is called get_ckn and located in the file hermite_expansion.py.

import sympy as sp
from sympy.printing.numpy import NumPyPrinter, \
    _known_functions_numpy, _known_constants_numpy
import os
from hermite_expansion import get_ckn

Afterwards, we define some symbols for SymPy

# Initialisation of symbolic variables
alpha, beta = sp.symbols('alpha beta', real=True, positive=True)
AX, BX = sp.symbols('A_x B_x', real=True)

as well as the overlap $E_{00}^{A_x, B_x}$

# Overlap for l_a = l_b = 0
S_00 = sp.sqrt(sp.pi / (alpha + beta)) * sp.exp(
    -((alpha * beta) / (alpha + beta)) * (AX**2 - 2 * AX * BX + BX**2)
)

Since we need Hermite Gaussian overlaps upto a certain maximum value of angular momentum, we shall write a function to generate them:

def generate_hermite_overlaps(lmax):
    hermite_overlaps = {}
    for k in range(0, lmax + 1):
        for l in range(0, lmax + 1):
            ho_kl = sp.simplify(sp.diff(sp.diff(S_00, AX, k), BX, l))
            hermite_overlaps[(k, l)] = ho_kl
    
    return hermite_overlaps

One specific overlap integral between two Cartesian Gaussians with angular momenta i and j can them be calculated using the following function:

def get_single_overlap(i, j, hermite_overlaps):
    overlap = 0
    for k in range(0, i + 1):
        cki = get_ckn(k, i, alpha)
        for l in range(0, j + 1):
            clj = get_ckn(l, j, beta)
            overlap += cki * clj * hermite_overlaps[(k, l)]
    overlap = sp.factor_terms(overlap)
    
    return overlap

We can then write a function to generate all Cartesian Gaussian overlaps upto a certain maximum angular momentum:

def generate_overlaps(lmax):
    hermite_overlaps = generate_hermite_overlaps(lmax)

    overlaps = {}
    # Loop through all combinations of Gaussian functions up to order lmax
    for i in range(lmax + 1):
        for j in range(lmax + 1):
            print(i, j)
            # Store the overlap integral in the dictionary with the key (i, j)
            overlaps[(i, j)] = get_single_overlap(i, j, hermite_overlaps)

    # Return the dictionary containing the overlap integrals
    return overlaps

We can then set

LMAX = 2

and generate formulas for all possible overlaps upto LMAX.

By inspecting the generated expressions, one might realize that some expressions occur very frequently, e.g. AX - BX, or alpha + beta. We can calculate these expressions once and store the value to avoid repeated calculations. To achieve this, we can substite these expressions with some new symbols:

# Substitute repeated expressions
X_AB, X_AB_SQ, P, Q = sp.symbols(
    'ab_diff ab_diff_squared ab_sum ab_product', 
    real=True,
)
subsdict = {
    AX - BX: X_AB,         
    AX**2 - 2 * AX * BX + BX**2: X_AB_SQ,
    alpha + beta: P, 
    alpha * beta: Q,
}

s_ij = generate_overlaps(LMAX)
s_ij = {k: v.subs(subsdict) for (k, v) in s_ij.items()}

We are almost done! The last step is to wrap these expressions into a Python function stored in a .py file. For this, we can write the following function:

def write_overlaps_py(overlaps, printer, path=''):
    with open(os.path.join(path, 'S.py'), 'w') as f:
        f.write('import numpy as np\n')
        f.write('def s_ij(i, j, alpha, beta, ax, bx):\n')
        # Calculate repeated expressions
        f.write('    ab_diff = ax - bx\n')
        f.write('    ab_diff_squared = ab_diff**2\n')
        f.write('    ab_sum = alpha + beta\n')
        f.write('    ab_product = alpha * beta\n')
        f.write('\n')
        # Write integrals for different cases
        
        for i, (key, value) in enumerate(overlaps.items()):
            if i == 0:
                if_str = 'if'
            else:
                if_str = 'elif'
            
            ia, ib = key
            code = printer.doprint(value)
            f.write(f'    {if_str} (i, j) == ({ia}, {ib}):\n')
            f.write(f'        return {code}\n')
        f.write('    else:\n')
        f.write('        raise NotImplementedError\n')

In this function, we have imported NumPy with the alias np. To convert the symbolic expressions into Python code with functions beginning with this alias, we setup a NumPyPrinter:

_numpy_known_functions = {k: f'np.{v}' for k, v 
                          in _known_functions_numpy.items()}
_numpy_known_constants = {k: f'np.{v}' for k, v 
                          in _known_constants_numpy.items()}

printer = NumPyPrinter()
printer._module = 'np'
printer.known_functions = _numpy_known_functions
printer.known_constants = _numpy_known_constants

And finally, we can generate the python file with all the integral expressions:

MY_PATH = '.'
write_overlaps_py(s_ij, printer, path=MY_PATH)

The path '.' stands for the location where you execute your python script. You can replace this with any valid path in your computer to generate S.py there.

Programming in Theoretical Chemistry

Overlap Integrals

Hermite Gaussian Overlaps

Cartesian Gaussian Overlaps

Code Generation