{ "cells": [ { "cell_type": "markdown", "id": "3d63629e", "metadata": {}, "source": [ "# Molecular Integrals in MiniHF\n", "\n", "This notebook provides a comprehensive introduction to molecular integrals in MiniHF. We'll cover:\n", "\n", "1. **Gaussian Type Orbitals (GTOs)** - The mathematical basis for quantum chemistry calculations\n", "2. **Basis Sets** - How GTOs are combined to form atomic basis functions\n", "3. **Obara-Saika Scheme** - The algorithm for computing molecular integrals efficiently\n", "4. **Practical Examples** - Hands-on demonstration of integral calculations in MiniHF\n", "\n", "Let's start by exploring the fundamental concepts." ] }, { "cell_type": "markdown", "id": "a0b7dfcb", "metadata": {}, "source": [ "## 1. Gaussian Type Orbitals (GTOs)\n", "\n", "MiniHF uses **Cartesian Gaussian Type Orbitals (GTOs)** as the basis functions for molecular calculations. A Cartesian GTO is defined mathematically as:\n", "\n", "$$\n", "G_{ijk}(\\mathbf{r}_A, \\alpha) = x_A^i y_A^j z_A^k e^{-\\alpha r_A^2}\n", "$$\n", "\n", "Where:\n", "- $\\mathbf{r}_A = \\mathbf{r} - \\mathbf{A}$ is the position vector relative to atom center $\\mathbf{A}$\n", "- $\\alpha > 0$ is the orbital exponent (controls the spatial extent)\n", "- $i, j, k \\geq 0$ are non-negative integers representing the angular momentum components\n", "- $l = i + j + k$ is the total angular momentum quantum number (0 = s, 1 = p, 2 = d, etc.)\n", "\n", "### Key Properties of GTOs\n", "\n", "1. **Factorization**: Cartesian GTOs factorize into three one-dimensional Gaussians:\n", " $$\n", " G_{ijk}(\\mathbf{r}_A, \\alpha) = G_i(x_A, \\alpha) G_j(y_A, \\alpha) G_k(z_A, \\alpha)\n", " $$\n", " where $G_i(x_A, \\alpha) = x_A^i e^{-\\alpha x_A^2}$\n", "\n", "2. **Gaussian Product Rule**: The product of two Gaussians is another Gaussian:\n", " $$\n", " e^{-\\alpha x_A^2} e^{-\\beta x_B^2} = e^{-\\mu X_{AB}^2} e^{-p x_P^2}\n", " $$\n", " where:\n", " - $p = \\alpha + \\beta$ (sum of exponents)\n", " - $\\mu = \\frac{\\alpha\\beta}{\\alpha+\\beta}$ (reduced exponent)\n", " - $P_x = \\frac{\\alpha A_x + \\beta B_x}{\\alpha+\\beta}$ (center-of-charge)\n", " - $X_{AB} = A_x - B_x$ (separation)\n", "\n", "3. **Differentiation**: The derivative of a Gaussian yields a linear combination of two Gaussians:\n", " $$\n", " \\frac{\\partial G_i(\\alpha, x_A)}{\\partial x} = 2\\alpha G_{i+1}(\\alpha, x_A) - i G_{i-1}(\\alpha, x_A)\n", " $$" ] }, { "cell_type": "markdown", "id": "basis-sets", "metadata": {}, "source": [ "## 2. Basis Sets\n", "\n", "Basis sets are collections of GTOs used to approximate atomic orbitals. Most quantum chemistry packages, including MiniHF, use the NWChem basis set format.\n", "\n", "### Basis Set Example: STO-3G\n", "\n", "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:\n", "\n", "```\n", "BASIS \"ao basis\" PRINT\n", "# Hydrogen: (3s) -> [1s]\n", "H S\n", " 3.42525091 0.15432897 \n", " 0.62391373 0.53532814 \n", " 0.16885540 0.44463454 \n", "# Oxygen: (6s,3p) -> [2s,1p]\n", "O S\n", " 130.7093200 0.15432897 \n", " 23.8088610 0.53532814 \n", " 6.4436083 0.44463454 \n", "O SP\n", " 5.0331513 -0.09996723 0.15591627 \n", " 1.1695961 0.39951283 0.60768372 \n", " 0.3803890 0.70011547 0.39195739 \n", "END\n", "```\n", "\n", "### Contracted Gaussian Type Orbitals (CGTOs)\n", "\n", "Contracted Gaussians are linear combinations of primitive Gaussians:\n", "$$\n", "\\chi_{\\mu}(\\mathbf{r}) = \\sum_{p=1}^P d_{\\mu p} G_{i,j,k}(\\mathbf{r}, \\alpha_p)\n", "$$\n", "where $d_{\\mu p}$ are contraction coefficients and $P$ is the number of primitives.\n", "\n", "### Available Basis Sets\n", "\n", "MiniHF supports basis sets from the EMSL Basis Set Exchange (https://bse.pnl.gov). Common basis sets include:\n", "- **STO-3G**: Minimal basis set\n", "- **6-31G**: Split-valence basis set\n", "- **6-311G**: Larger split-valence basis set\n", "- **cc-pVDZ/cc-pVTZ**: Correlation-consistent basis sets" ] }, { "cell_type": "markdown", "id": "obara-saika", "metadata": {}, "source": [ "## 3. The Obara-Saika Recurrence Scheme\n", "\n", "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.\n", "\n", "### 3.1 Overlap Integrals\n", "\n", "The overlap integral measures the overlap between two basis functions:\n", "$$\n", "S_{ab} = \\langle G_a | G_b \\rangle\n", "$$\n", "\n", "It factorizes into three Cartesian components:\n", "$$\n", "S_{ab} = S_{ij}(x) S_{kl}(y) S_{mn}(z)\n", "$$\n", "\n", "**Recurrence Relation**:\n", "$$\n", "\\begin{aligned}\n", "S_{i+1,j} &= X_{PA} S_{ij} + \\frac{1}{2p}(i S_{i-1,j} + j S_{i,j-1}) \\\\\n", "S_{i,j+1} &= X_{PB} S_{ij} + \\frac{1}{2p}(i S_{i-1,j} + j S_{i,j-1})\n", "\\end{aligned}\n", "$$\n", "\n", "**Boundary Condition**:\n", "$$\n", "S_{0,0} = \\sqrt{\\frac{\\pi}{p}} e^{-\\mu X_{AB}^2}\n", "$$" ] }, { "cell_type": "markdown", "id": "kinetic-integrals", "metadata": {}, "source": [ "### 3.2 Kinetic Energy Integrals\n", "\n", "The kinetic energy integral describes the electron's kinetic energy:\n", "$$\n", "T_{ab} = \\langle G_a | -\\frac{1}{2} \\nabla^2 | G_b \\rangle\n", "$$\n", "\n", "It also factorizes into Cartesian components:\n", "$$\n", "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)\n", "$$\n", "\n", "**Recurrence Relation**:\n", "$$\n", "\\begin{aligned}\n", "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}) \\\\\n", "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})\n", "\\end{aligned}\n", "$$\n", "\n", "**Boundary Condition**:\n", "$$\n", "T_{0,0} = \\left(a - 2a^2 \\left(X_{PA}^2 + \\frac{1}{2p}\\right)\\right) S_{0,0}\n", "$$" ] }, { "cell_type": "markdown", "id": "nuclear-attraction", "metadata": {}, "source": [ "### 3.3 Nuclear Attraction Integrals\n", "\n", "The nuclear attraction integral describes the electrostatic attraction between electrons and nuclei:\n", "$$\n", "V_{ab} = \\langle G_a | -\\sum_{\\alpha} \\frac{Z_{\\alpha}}{|\\mathbf{r} - \\mathbf{R}_{\\alpha}|} | G_b \\rangle\n", "$$\n", "\n", "**Recurrence Relation** (Obara-Saika):\n", "$$\n", "\\begin{aligned}\n", "\\Theta_{i+1,j,k,l,m,n}^{N} &= X_{PA} \\Theta_{i,j,k,l,m,n}^{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}) \\\\\n", "&- X_{PC} \\Theta_{i,j,k,l,m,n}^{N+1} \\\\\n", "&- \\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})\n", "\\end{aligned}\n", "$$\n", "\n", "**Boundary Condition**:\n", "$$\n", "\\Theta_{000000}^{N} = \\frac{2\\pi}{p} K_{ab}^{xyz} F_N(p R_{PC}^2)\n", "$$\n", "where $F_N(x)$ is the Boys function and $K_{ab}^{xyz} = e^{-\\mu R_{AB}^2}$." ] }, { "cell_type": "markdown", "id": "electron-repulsion", "metadata": {}, "source": [ "### 3.4 Electron Repulsion Integrals (ERIs)\n", "\n", "The electron repulsion integral describes the Coulomb repulsion between pairs of electrons:\n", "$$\n", "(ab|cd) = \\langle G_a G_b | \\frac{1}{|\\mathbf{r}_1 - \\mathbf{r}_2|} | G_c G_d \\rangle\n", "$$\n", "\n", "This is the most computationally expensive integral in Hartree-Fock theory, scaling as $O(N^4)$ with the number of basis functions.\n", "\n", "**Boundary Condition**:\n", "$$\n", "\\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)\n", "$$\n", "\n", "**Two-Electron Recurrence Relation**:\n", "$$\n", "\\begin{aligned}\n", "\\Theta_{i+1,j,k,l}^{N} &= X_{PA} \\Theta_{i,j,k,l}^{N} \\\\\n", "&- \\frac{\\alpha}{p} X_{PQ} \\Theta_{i,j,k,l}^{N+1} \\\\\n", "&+ \\frac{i}{2p} \\left(\\Theta_{i-1,j,k,l}^{N} - \\frac{\\alpha}{p} \\Theta_{i-1,j,k,l}^{N+1}\\right) \\\\\n", "&+ \\frac{j}{2p} \\left(\\Theta_{i,j-1,k,l}^{N} - \\frac{\\alpha}{p} \\Theta_{i,j-1,k,l}^{N+1}\\right) \\\\\n", "&- \\frac{k}{2(p+q)} \\Theta_{i,j,k-1,l}^{N+1} \\\\\n", "&- \\frac{l}{2(p+q)} \\Theta_{i,j,k,l-1}^{N+1}\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "id": "practical-examples", "metadata": {}, "source": [ "## 4. Practical Examples with MiniHF\n", "\n", "Let's see how to compute molecular integrals using MiniHF. We'll create a simple molecule and calculate the Hamiltonian matrix." ] }, { "cell_type": "code", "execution_count": 1, "id": "create-molecule", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Molecule: Molecule(atoms=2, charge=0, multiplicity=1)\n", "Number of atoms: 2\n", "Number of basis functions: 2\n", "Basis function types:\n", " 0: (0, 0, 0) (l=0)\n", " 1: (0, 0, 0) (l=0)\n" ] } ], "source": [ "# Import required modules\n", "from minihf.molecule import Molecule\n", "from minihf.basis import BasisSet\n", "from minihf.hamiltonian import Hamiltonian\n", "import numpy as np\n", "\n", "# Create a simple hydrogen molecule\n", "atoms = ['H 0.0, 0.0, -0.37', 'H 0.0, 0.0, 0.37']\n", "mol = Molecule.build(atoms=atoms, charge=0, multiplicity=1)\n", "\n", "# Build basis set (STO-3G for minimal example)\n", "basis = BasisSet.build(mol, 'sto-3g.nwchem')\n", "\n", "print(f\"Molecule: {mol}\")\n", "print(f\"Number of atoms: {len(mol)}\")\n", "print(f\"Number of basis functions: {len(basis)}\")\n", "print(f\"Basis function types:\")\n", "for i, orb in enumerate(basis):\n", " print(f\" {i}: {orb.shell} (l={sum(orb.shell)})\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "compute-integrals", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Integral Matrix Dimensions ===\n", "Overlap matrix (S): (2, 2)\n", "Kinetic matrix (T): (2, 2)\n", "Potential matrix (V): (2, 2)\n", "Core Hamiltonian (Hcore): (2, 2)\n", "Electron repulsion integrals (Eri): (2, 2, 2, 2)\n" ] } ], "source": [ "# Compute the Hamiltonian\n", "ham = Hamiltonian.build(mol, basis)\n", "\n", "# Extract integral matrices\n", "S = ham['S'] # Overlap matrix\n", "T = ham['T'] # Kinetic energy matrix\n", "V = ham['V'] # Potential energy matrix\n", "Hcore = ham['Hcore'] # Core Hamiltonian (T + V)\n", "Eri = ham['Eri'] # Electron repulsion integrals\n", "\n", "# Print matrix dimensions\n", "print(\"=== Integral Matrix Dimensions ===\")\n", "print(f\"Overlap matrix (S): {S.shape}\")\n", "print(f\"Kinetic matrix (T): {T.shape}\")\n", "print(f\"Potential matrix (V): {V.shape}\")\n", "print(f\"Core Hamiltonian (Hcore): {Hcore.shape}\")\n", "print(f\"Electron repulsion integrals (Eri): {Eri.shape}\")" ] }, { "cell_type": "code", "execution_count": 3, "id": "analyze-overlap", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Overlap Matrix Analysis ===\n", "\n", "Overlap matrix S:\n", "[[1. 0.878182]\n", " [0.878182 1. ]]\n", "\n", "Trace of S: 2.000000\n", "Determinant of S: 0.228797\n", "S is symmetric: True\n", "\n", "Eigenvalues of S: [0.121818 1.878182]\n" ] } ], "source": [ "# Analyze the overlap matrix\n", "print(\"=== Overlap Matrix Analysis ===\")\n", "print(f\"\\nOverlap matrix S:\")\n", "print(np.round(S, 6))\n", "\n", "# Check key properties\n", "print(f\"\\nTrace of S: {np.trace(S):.6f}\")\n", "print(f\"Determinant of S: {np.linalg.det(S):.6f}\")\n", "print(f\"S is symmetric: {np.allclose(S, S.T)}\")\n", "\n", "# Diagonalize S to check eigenvalues\n", "eigvals = np.linalg.eigvalsh(S)\n", "print(f\"\\nEigenvalues of S: {np.round(eigvals, 6)}\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "analyze-core-hamiltonian", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Core Hamiltonian Analysis ===\n", "\n", "Core Hamiltonian Hcore = T + V:\n", "[[-1.404011 -1.391124]\n", " [-1.391124 -1.404011]]\n", "\n", "Verification: |Hcore - (T + V)|_max = 0.000000e+00\n", "Hcore is symmetric: True\n", "\n", "Diagonal elements of Hcore:\n", " Hcore[0,0] = -1.404011\n", " Hcore[1,1] = -1.404011\n" ] } ], "source": [ "# Analyze the core Hamiltonian\n", "print(\"=== Core Hamiltonian Analysis ===\")\n", "print(f\"\\nCore Hamiltonian Hcore = T + V:\")\n", "print(np.round(Hcore, 6))\n", "\n", "# Verify that Hcore = T + V\n", "print(f\"\\nVerification: |Hcore - (T + V)|_max = {np.max(np.abs(Hcore - (T + V))):.6e}\")\n", "\n", "# Check symmetry\n", "print(f\"Hcore is symmetric: {np.allclose(Hcore, Hcore.T)}\")\n", "\n", "# Get diagonal elements (approximate orbital energies)\n", "print(f\"\\nDiagonal elements of Hcore:\")\n", "for i in range(len(basis)):\n", " print(f\" Hcore[{i},{i}] = {Hcore[i,i]:.6f}\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "analyze-eri", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Electron Repulsion Integrals Analysis ===\n", "Total ERIs: 16\n", "Unique ERIs (with symmetry): 9\n", "\n", "Selected ERI values:\n", " (0,0|0,0) = 0.774606\n", " (0,0|1,1) = 0.699480\n", " (0,1|0,1) = 0.570752\n", " (0,1|1,0) = 0.570752\n", "\n", "ERI Symmetry Checks:\n", " (0,0|1,1) == (1,1|0,0): True\n", " (0,1|0,1) == (1,0|1,0): True\n", " (0,1|0,1) == (0,1|1,0): True\n" ] } ], "source": [ "# Analyze electron repulsion integrals\n", "print(\"=== Electron Repulsion Integrals Analysis ===\")\n", "\n", "# Number of unique ERIs (8-fold symmetry)\n", "n_bf = len(basis)\n", "total_eris = n_bf ** 4\n", "unique_eris = (n_bf * (n_bf + 1) // 2) ** 2\n", "print(f\"Total ERIs: {total_eris}\")\n", "print(f\"Unique ERIs (with symmetry): {unique_eris}\")\n", "\n", "# Check some specific ERIs\n", "print(\"\\nSelected ERI values:\")\n", "print(f\" (0,0|0,0) = {Eri[0,0,0,0]:.6f}\")\n", "print(f\" (0,0|1,1) = {Eri[0,0,1,1]:.6f}\")\n", "print(f\" (0,1|0,1) = {Eri[0,1,0,1]:.6f}\")\n", "print(f\" (0,1|1,0) = {Eri[0,1,1,0]:.6f}\")\n", "\n", "# Verify symmetry relations\n", "print(\"\\nERI Symmetry Checks:\")\n", "print(f\" (0,0|1,1) == (1,1|0,0): {np.isclose(Eri[0,0,1,1], Eri[1,1,0,0])}\")\n", "print(f\" (0,1|0,1) == (1,0|1,0): {np.isclose(Eri[0,1,0,1], Eri[1,0,1,0])}\")\n", "print(f\" (0,1|0,1) == (0,1|1,0): {np.isclose(Eri[0,1,0,1], Eri[0,1,1,0])}\")" ] }, { "cell_type": "markdown", "id": "summary", "metadata": {}, "source": [ "## 5. Summary\n", "\n", "In this notebook, we've covered:\n", "\n", "1. **Gaussian Type Orbitals**: The mathematical foundation of modern quantum chemistry\n", "2. **Basis Sets**: Collections of GTOs used to approximate atomic orbitals\n", "3. **Obara-Saika Algorithm**: Efficient recursive computation of molecular integrals\n", "4. **Practical Usage**: How to compute integrals in MiniHF\n", "\n", "### Key Takeaways\n", "\n", "- **Overlap integrals** measure how much two basis functions overlap\n", "- **Kinetic energy integrals** describe electron motion\n", "- **Nuclear attraction integrals** describe electron-nucleus interactions\n", "- **Electron repulsion integrals** describe electron-electron interactions (most expensive!)\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "id": "bonus", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Water Molecule Integrals ===\n", "Water molecule with 3 atoms\n", "Basis functions: 13\n", "ERI tensor size: (13, 13, 13, 13)\n", "Total ERIs: 28,561\n", "\n", "Overlap matrix trace: 13.000000\n", "Core Hamiltonian trace: -123.119025\n", "Smallest Hcore eigenvalue: -48.772084\n", "Largest Hcore eigenvalue: -0.075431\n" ] } ], "source": [ "# Bonus: Compute integrals for water molecule\n", "print(\"=== Water Molecule Integrals ===\")\n", "\n", "# Create water molecule\n", "h2o_atoms = ['O 0.0, 0.0, 0.0', 'H 0.0, -0.757, 0.587', 'H 0.0, 0.757, 0.587']\n", "h2o = Molecule.build(atoms=h2o_atoms, charge=0, multiplicity=1)\n", "\n", "# Build 6-31G basis set\n", "h2o_basis = BasisSet.build(h2o, '6-31g.nwchem')\n", "\n", "# Compute Hamiltonian\n", "h2o_ham = Hamiltonian.build(h2o, h2o_basis)\n", "\n", "print(f\"Water molecule with {len(h2o)} atoms\")\n", "print(f\"Basis functions: {len(h2o_basis)}\")\n", "print(f\"ERI tensor size: {h2o_ham['Eri'].shape}\")\n", "print(f\"Total ERIs: {h2o_ham['Eri'].size:,}\")\n", "\n", "# Compute some properties\n", "h2o_S = h2o_ham['S']\n", "h2o_Hcore = h2o_ham['Hcore']\n", "\n", "print(f\"\\nOverlap matrix trace: {np.trace(h2o_S):.6f}\")\n", "print(f\"Core Hamiltonian trace: {np.trace(h2o_Hcore):.6f}\")\n", "print(f\"Smallest Hcore eigenvalue: {np.linalg.eigvalsh(h2o_Hcore).min():.6f}\")\n", "print(f\"Largest Hcore eigenvalue: {np.linalg.eigvalsh(h2o_Hcore).max():.6f}\")" ] } ], "metadata": { "kernelspec": { "display_name": "base", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.12" } }, "nbformat": 4, "nbformat_minor": 5 }