Molecular Integrals in MiniHF

This notebook provides a comprehensive introduction to molecular integrals in MiniHF. We’ll cover:

  1. Gaussian Type Orbitals (GTOs) - The mathematical basis for quantum chemistry calculations

  2. Basis Sets - How GTOs are combined to form atomic basis functions

  3. Obara-Saika Scheme - The algorithm for computing molecular integrals efficiently

  4. Practical Examples - Hands-on demonstration of integral calculations in MiniHF

Let’s start by exploring the fundamental concepts.

1. Gaussian Type Orbitals (GTOs)

MiniHF uses Cartesian Gaussian Type Orbitals (GTOs) as the basis functions for molecular calculations. A Cartesian GTO is defined mathematically as:

\[ G_{ijk}(\mathbf{r}_A, \alpha) = x_A^i y_A^j z_A^k e^{-\alpha r_A^2} \]

Where:

  • \(\mathbf{r}_A = \mathbf{r} - \mathbf{A}\) is the position vector relative to atom center \(\mathbf{A}\)

  • \(\alpha > 0\) is the orbital exponent (controls the spatial extent)

  • \(i, j, k \geq 0\) are non-negative integers representing the angular momentum components

  • \(l = i + j + k\) is the total angular momentum quantum number (0 = s, 1 = p, 2 = d, etc.)

Key Properties of GTOs

  1. Factorization: Cartesian GTOs factorize into three one-dimensional Gaussians: $\( G_{ijk}(\mathbf{r}_A, \alpha) = G_i(x_A, \alpha) G_j(y_A, \alpha) G_k(z_A, \alpha) \)\( where \)G_i(x_A, \alpha) = x_A^i e^{-\alpha x_A^2}$

  2. Gaussian Product Rule: The product of two Gaussians is another Gaussian: $\( e^{-\alpha x_A^2} e^{-\beta x_B^2} = e^{-\mu X_{AB}^2} e^{-p x_P^2} \)$ where:

    • \(p = \alpha + \beta\) (sum of exponents)

    • \(\mu = \frac{\alpha\beta}{\alpha+\beta}\) (reduced exponent)

    • \(P_x = \frac{\alpha A_x + \beta B_x}{\alpha+\beta}\) (center-of-charge)

    • \(X_{AB} = A_x - B_x\) (separation)

  3. Differentiation: The derivative of a Gaussian yields a linear combination of two Gaussians: $\( \frac{\partial G_i(\alpha, x_A)}{\partial x} = 2\alpha G_{i+1}(\alpha, x_A) - i G_{i-1}(\alpha, x_A) \)$

2. Basis Sets

Basis sets are collections of GTOs used to approximate atomic orbitals. Most quantum chemistry packages, including MiniHF, use the NWChem basis set format.

Basis Set Example: STO-3G

The STO-3G basis set approximates Slater-type orbitals (STOs) with 3 primitive Gaussians (PGTOs) contracted together. Here’s the definition for hydrogen and oxygen:

BASIS "ao basis" PRINT
# Hydrogen: (3s) -> [1s]
H    S
    3.42525091             0.15432897       
    0.62391373             0.53532814       
    0.16885540             0.44463454       
# Oxygen: (6s,3p) -> [2s,1p]
O    S
    130.7093200              0.15432897       
    23.8088610              0.53532814       
    6.4436083              0.44463454       
O    SP
    5.0331513             -0.09996723             0.15591627       
    1.1695961              0.39951283             0.60768372       
    0.3803890              0.70011547             0.39195739       
END

Contracted Gaussian Type Orbitals (CGTOs)

Contracted Gaussians are linear combinations of primitive Gaussians: $\( \chi_{\mu}(\mathbf{r}) = \sum_{p=1}^P d_{\mu p} G_{i,j,k}(\mathbf{r}, \alpha_p) \)\( where \)d_{\mu p}\( are contraction coefficients and \)P$ is the number of primitives.

Available Basis Sets

MiniHF supports basis sets from the EMSL Basis Set Exchange (https://bse.pnl.gov). Common basis sets include:

  • STO-3G: Minimal basis set

  • 6-31G: Split-valence basis set

  • 6-311G: Larger split-valence basis set

  • cc-pVDZ/cc-pVTZ: Correlation-consistent basis sets

3. The Obara-Saika Recurrence Scheme

The Obara-Saika algorithm is a recursive method for efficiently computing molecular integrals over Gaussian basis functions. It avoids direct evaluation of complicated integrals by using recurrence relations.

3.1 Overlap Integrals

The overlap integral measures the overlap between two basis functions: $\( S_{ab} = \langle G_a | G_b \rangle \)$

It factorizes into three Cartesian components: $\( S_{ab} = S_{ij}(x) S_{kl}(y) S_{mn}(z) \)$

Recurrence Relation: $\( \begin{aligned} S_{i+1,j} &= X_{PA} S_{ij} + \frac{1}{2p}(i S_{i-1,j} + j S_{i,j-1}) \\ S_{i,j+1} &= X_{PB} S_{ij} + \frac{1}{2p}(i S_{i-1,j} + j S_{i,j-1}) \end{aligned} \)$

Boundary Condition: $\( S_{0,0} = \sqrt{\frac{\pi}{p}} e^{-\mu X_{AB}^2} \)$

3.2 Kinetic Energy Integrals

The kinetic energy integral describes the electron’s kinetic energy: $\( T_{ab} = \langle G_a | -\frac{1}{2} \nabla^2 | G_b \rangle \)$

It also factorizes into Cartesian components: $\( T_{ab} = T_{ij}(x) S_{kl}(y) S_{mn}(z) + S_{ij}(x) T_{kl}(y) S_{mn}(z) + S_{ij}(x) S_{kl}(y) T_{mn}(z) \)$

Recurrence Relation: $\( \begin{aligned} T_{i+1,j} &= X_{PA} T_{ij} + \frac{1}{2p}(i T_{i-1,j} + j T_{i,j-1}) + \frac{b}{p}(2a S_{i+1,j} - i S_{i-1,j}) \\ T_{i,j+1} &= X_{PB} T_{ij} + \frac{1}{2p}(i T_{i-1,j} + j T_{i,j-1}) + \frac{a}{p}(2b S_{i,j+1} - j S_{i,j-1}) \end{aligned} \)$

Boundary Condition: $\( T_{0,0} = \left(a - 2a^2 \left(X_{PA}^2 + \frac{1}{2p}\right)\right) S_{0,0} \)$

3.3 Nuclear Attraction Integrals

The nuclear attraction integral describes the electrostatic attraction between electrons and nuclei: $\( V_{ab} = \langle G_a | -\sum_{\alpha} \frac{Z_{\alpha}}{|\mathbf{r} - \mathbf{R}_{\alpha}|} | G_b \rangle \)$

Recurrence Relation (Obara-Saika): $\( \begin{aligned} \Theta_{i+1,j,k,l,m,n}^{N} &= X_{PA} \Theta_{i,j,k,l,m,n}^{N} \\ &+ \frac{1}{2p}(i \Theta_{i-1,j,k,l,m,n}^{N} + j \Theta_{i,j-1,k,l,m,n}^{N}) \\ &- X_{PC} \Theta_{i,j,k,l,m,n}^{N+1} \\ &- \frac{1}{2p}(i \Theta_{i-1,j,k,l,m,n}^{N+1} + j \Theta_{i,j-1,k,l,m,n}^{N+1}) \end{aligned} \)$

Boundary Condition: $\( \Theta_{000000}^{N} = \frac{2\pi}{p} K_{ab}^{xyz} F_N(p R_{PC}^2) \)\( where \)F_N(x)\( is the Boys function and \)K_{ab}^{xyz} = e^{-\mu R_{AB}^2}$.

3.4 Electron Repulsion Integrals (ERIs)

The electron repulsion integral describes the Coulomb repulsion between pairs of electrons: $\( (ab|cd) = \langle G_a G_b | \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} | G_c G_d \rangle \)$

This is the most computationally expensive integral in Hartree-Fock theory, scaling as \(O(N^4)\) with the number of basis functions.

Boundary Condition: $\( \Theta_{0000;0000;0000}^{N} = \frac{2\pi^{2.5}}{pq\sqrt{p+q}} K_{ab}^{xyz} K_{cd}^{xyz} F_N(\alpha R_{PQ}^2) \)$

Two-Electron Recurrence Relation: $\( \begin{aligned} \Theta_{i+1,j,k,l}^{N} &= X_{PA} \Theta_{i,j,k,l}^{N} \\ &- \frac{\alpha}{p} X_{PQ} \Theta_{i,j,k,l}^{N+1} \\ &+ \frac{i}{2p} \left(\Theta_{i-1,j,k,l}^{N} - \frac{\alpha}{p} \Theta_{i-1,j,k,l}^{N+1}\right) \\ &+ \frac{j}{2p} \left(\Theta_{i,j-1,k,l}^{N} - \frac{\alpha}{p} \Theta_{i,j-1,k,l}^{N+1}\right) \\ &- \frac{k}{2(p+q)} \Theta_{i,j,k-1,l}^{N+1} \\ &- \frac{l}{2(p+q)} \Theta_{i,j,k,l-1}^{N+1} \end{aligned} \)$

4. Practical Examples with MiniHF

Let’s see how to compute molecular integrals using MiniHF. We’ll create a simple molecule and calculate the Hamiltonian matrix.

# Import required modules
from minihf.molecule import Molecule
from minihf.basis import BasisSet
from minihf.hamiltonian import Hamiltonian
import numpy as np

# Create a simple hydrogen molecule
atoms = ['H 0.0, 0.0, -0.37', 'H 0.0, 0.0, 0.37']
mol = Molecule.build(atoms=atoms, charge=0, multiplicity=1)

# Build basis set (STO-3G for minimal example)
basis = BasisSet.build(mol, 'sto-3g.nwchem')

print(f"Molecule: {mol}")
print(f"Number of atoms: {len(mol)}")
print(f"Number of basis functions: {len(basis)}")
print(f"Basis function types:")
for i, orb in enumerate(basis):
    print(f"  {i}: {orb.shell} (l={sum(orb.shell)})")

Molecule: Molecule(atoms=2, charge=0, multiplicity=1)
Number of atoms: 2
Number of basis functions: 2
Basis function types:
  0: (0, 0, 0) (l=0)
  1: (0, 0, 0) (l=0)

# Compute the Hamiltonian
ham = Hamiltonian.build(mol, basis)

# Extract integral matrices
S = ham['S']  # Overlap matrix
T = ham['T']  # Kinetic energy matrix
V = ham['V']  # Potential energy matrix
Hcore = ham['Hcore']  # Core Hamiltonian (T + V)
Eri = ham['Eri']  # Electron repulsion integrals

# Print matrix dimensions
print("=== Integral Matrix Dimensions ===")
print(f"Overlap matrix (S): {S.shape}")
print(f"Kinetic matrix (T): {T.shape}")
print(f"Potential matrix (V): {V.shape}")
print(f"Core Hamiltonian (Hcore): {Hcore.shape}")
print(f"Electron repulsion integrals (Eri): {Eri.shape}")

=== Integral Matrix Dimensions ===
Overlap matrix (S): (2, 2)
Kinetic matrix (T): (2, 2)
Potential matrix (V): (2, 2)
Core Hamiltonian (Hcore): (2, 2)
Electron repulsion integrals (Eri): (2, 2, 2, 2)

# Analyze the overlap matrix
print("=== Overlap Matrix Analysis ===")
print(f"\nOverlap matrix S:")
print(np.round(S, 6))

# Check key properties
print(f"\nTrace of S: {np.trace(S):.6f}")
print(f"Determinant of S: {np.linalg.det(S):.6f}")
print(f"S is symmetric: {np.allclose(S, S.T)}")

# Diagonalize S to check eigenvalues
eigvals = np.linalg.eigvalsh(S)
print(f"\nEigenvalues of S: {np.round(eigvals, 6)}")

=== Overlap Matrix Analysis ===

Overlap matrix S:
[[1.       0.878182]
 [0.878182 1.      ]]

Trace of S: 2.000000
Determinant of S: 0.228797
S is symmetric: True

Eigenvalues of S: [0.121818 1.878182]

# Analyze the core Hamiltonian
print("=== Core Hamiltonian Analysis ===")
print(f"\nCore Hamiltonian Hcore = T + V:")
print(np.round(Hcore, 6))

# Verify that Hcore = T + V
print(f"\nVerification: |Hcore - (T + V)|_max = {np.max(np.abs(Hcore - (T + V))):.6e}")

# Check symmetry
print(f"Hcore is symmetric: {np.allclose(Hcore, Hcore.T)}")

# Get diagonal elements (approximate orbital energies)
print(f"\nDiagonal elements of Hcore:")
for i in range(len(basis)):
    print(f"  Hcore[{i},{i}] = {Hcore[i,i]:.6f}")

=== Core Hamiltonian Analysis ===

Core Hamiltonian Hcore = T + V:
[[-1.404011 -1.391124]
 [-1.391124 -1.404011]]

Verification: |Hcore - (T + V)|_max = 0.000000e+00
Hcore is symmetric: True

Diagonal elements of Hcore:
  Hcore[0,0] = -1.404011
  Hcore[1,1] = -1.404011

# Analyze electron repulsion integrals
print("=== Electron Repulsion Integrals Analysis ===")

# Number of unique ERIs (8-fold symmetry)
n_bf = len(basis)
total_eris = n_bf ** 4
unique_eris = (n_bf * (n_bf + 1) // 2) ** 2
print(f"Total ERIs: {total_eris}")
print(f"Unique ERIs (with symmetry): {unique_eris}")

# Check some specific ERIs
print("\nSelected ERI values:")
print(f"  (0,0|0,0) = {Eri[0,0,0,0]:.6f}")
print(f"  (0,0|1,1) = {Eri[0,0,1,1]:.6f}")
print(f"  (0,1|0,1) = {Eri[0,1,0,1]:.6f}")
print(f"  (0,1|1,0) = {Eri[0,1,1,0]:.6f}")

# Verify symmetry relations
print("\nERI Symmetry Checks:")
print(f"  (0,0|1,1) == (1,1|0,0): {np.isclose(Eri[0,0,1,1], Eri[1,1,0,0])}")
print(f"  (0,1|0,1) == (1,0|1,0): {np.isclose(Eri[0,1,0,1], Eri[1,0,1,0])}")
print(f"  (0,1|0,1) == (0,1|1,0): {np.isclose(Eri[0,1,0,1], Eri[0,1,1,0])}")

=== Electron Repulsion Integrals Analysis ===
Total ERIs: 16
Unique ERIs (with symmetry): 9

Selected ERI values:
  (0,0|0,0) = 0.774606
  (0,0|1,1) = 0.699480
  (0,1|0,1) = 0.570752
  (0,1|1,0) = 0.570752

ERI Symmetry Checks:
  (0,0|1,1) == (1,1|0,0): True
  (0,1|0,1) == (1,0|1,0): True
  (0,1|0,1) == (0,1|1,0): True

5. Summary

In this notebook, we’ve covered:

  1. Gaussian Type Orbitals: The mathematical foundation of modern quantum chemistry

  2. Basis Sets: Collections of GTOs used to approximate atomic orbitals

  3. Obara-Saika Algorithm: Efficient recursive computation of molecular integrals

  4. Practical Usage: How to compute integrals in MiniHF

Key Takeaways

  • Overlap integrals measure how much two basis functions overlap

  • Kinetic energy integrals describe electron motion

  • Nuclear attraction integrals describe electron-nucleus interactions

  • Electron repulsion integrals describe electron-electron interactions (most expensive!)

These integrals form the building blocks for Hartree-Fock and post-Hartree-Fock methods. The Obara-Saika recurrence scheme makes it feasible to compute these integrals efficiently for large basis sets.

# Bonus: Compute integrals for water molecule
print("=== Water Molecule Integrals ===")

# Create water molecule
h2o_atoms = ['O 0.0, 0.0, 0.0', 'H 0.0, -0.757, 0.587', 'H 0.0, 0.757, 0.587']
h2o = Molecule.build(atoms=h2o_atoms, charge=0, multiplicity=1)

# Build 6-31G basis set
h2o_basis = BasisSet.build(h2o, '6-31g.nwchem')

# Compute Hamiltonian
h2o_ham = Hamiltonian.build(h2o, h2o_basis)

print(f"Water molecule with {len(h2o)} atoms")
print(f"Basis functions: {len(h2o_basis)}")
print(f"ERI tensor size: {h2o_ham['Eri'].shape}")
print(f"Total ERIs: {h2o_ham['Eri'].size:,}")

# Compute some properties
h2o_S = h2o_ham['S']
h2o_Hcore = h2o_ham['Hcore']

print(f"\nOverlap matrix trace: {np.trace(h2o_S):.6f}")
print(f"Core Hamiltonian trace: {np.trace(h2o_Hcore):.6f}")
print(f"Smallest Hcore eigenvalue: {np.linalg.eigvalsh(h2o_Hcore).min():.6f}")
print(f"Largest Hcore eigenvalue: {np.linalg.eigvalsh(h2o_Hcore).max():.6f}")

=== Water Molecule Integrals ===
Water molecule with 3 atoms
Basis functions: 13
ERI tensor size: (13, 13, 13, 13)
Total ERIs: 28,561

Overlap matrix trace: 13.000000
Core Hamiltonian trace: -123.119025
Smallest Hcore eigenvalue: -48.772084
Largest Hcore eigenvalue: -0.075431