H2 Atom and FCI

Full Configuration Interaction

The H2 molecule is a fundamental system for understanding Full Configuration Interaction (FCI) in quantum chemistry. With two electrons and a minimal basis of two 1s atomic orbitals (one on each hydrogen atom), H2 provides a simple yet instructive example of FCI that can be solved exactly.

System Setup

For H2, we have:

2 electrons (one from each hydrogen atom)
2 spatial orbitals: 1s orbitals centered on each H atom ( $ϕ_{a}$ and $ϕ_{b}$ )
4 spin orbitals: $ϕ_{a ↑}, ϕ_{a ↓}, ϕ_{b ↑}, ϕ_{b ↓}$

FCI Configuration Space

For 2 electrons in 4 spin orbitals, the number of configurations is:

N_{config} = (2 4) = 6

The six possible Slater determinants are:

$∣ Φ_{1} ⟩ = ∣ ϕ_{a ↑} ϕ_{a ↓} ⟩$ - Both electrons on atom A
$∣ Φ_{2} ⟩ = ∣ ϕ_{b ↑} ϕ_{b ↓} ⟩$ - Both electrons on atom B
$∣ Φ_{3} ⟩ = ∣ ϕ_{a ↑} ϕ_{b ↑} ⟩$ - Both electrons spin up, one on each atom
$∣ Φ_{4} ⟩ = ∣ ϕ_{a ↓} ϕ_{b ↓} ⟩$ - Both electrons spin down, one on each atom
$∣ Φ_{5} ⟩ = ∣ ϕ_{a ↑} ϕ_{b ↓} ⟩$ - One electron on each atom, opposite spins
$∣ Φ_{6} ⟩ = ∣ ϕ_{a ↓} ϕ_{b ↑} ⟩$ - One electron on each atom, opposite spins (different ordering)

Hamiltonian

The electronic Hamiltonian for H2 is:

\hat{H} = i = 1 \sum 2 \hat{h} (i) + \frac{1}{r _{12}} + \frac{1}{R}

where:

$\hat{h} (i) = - \frac{1}{2} \nabla_{i}^{2} - \frac{1}{r _{i A}} - \frac{1}{r _{i B}}$ is the one-electron operator
$\frac{1}{r _{12}}$ is the electron-electron repulsion
$\frac{1}{R}$ is the nuclear-nuclear repulsion (constant for fixed bond length)

Hamiltonian Matrix Elements

Using the Slater-Condon rules, the Hamiltonian matrix elements are computed as:

H_{I J} = ⟨ Φ_{I} ∣ \hat{H} ∣ Φ_{J} ⟩

For determinants that differ by 0, 1, or 2 orbitals, we use:

Same determinant: $H_{II} = \sum_{i} h_{ii} + \sum_{i < j} (ii ∣ jj) - (ij ∣ ji)$
Differ by one orbital: $H_{I J} = h_{pq} + \sum_{r} [(p r ∣ q r) - (p r ∣ r q)]$
Differ by two orbitals: $H_{I J} = (pq ∣ rs) - (pq ∣ sr)$

where $(pq ∣ rs)$ denotes the two-electron integral in chemist's notation.

Integral Evaluation

For 1s Gaussian-type orbitals, we need to compute:

Overlap: $S_{ab} = ⟨ ϕ_{a} ∣ ϕ_{b} ⟩$
Kinetic energy: $T_{ab} = ⟨ ϕ_{a} ∣ - \frac{1}{2} \nabla^{2} ∣ ϕ_{b} ⟩$
Electron-nucleus: $V_{ab}^{A} = ⟨ ϕ_{a} ∣ - \frac{1}{r _{A}} ∣ ϕ_{b} ⟩$ and $V_{ab}^{B}$
Two-electron: $(ab ∣ c d) = \int\int ϕ_{a} (1) ϕ_{b} (1) \frac{1}{r _{12}} ϕ_{c} (2) ϕ_{d} (2) d r_{1} d r_{2}$

These integrals can be computed analytically using Gaussian-type orbitals. See the 1s Integrals page for formulas and implementations.

Python Implementation

Here is a complete Python implementation to solve H2 using FCI:

import numpy as np
from scipy.linalg import eigh
from scipy.special import erf

def F0(t):
    """Compute the F0 function used in electron integrals."""
    if t < 1e-10:
        return 1.0
    sqrt_t = np.sqrt(t)
    return np.sqrt(np.pi / (4 * t)) * erf(sqrt_t)

def overlap(alpha, RA, beta, RB):
    """Overlap integral for 1s Gaussians."""
    RA = np.array(RA)
    RB = np.array(RB)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    return (np.pi / p) ** (3.0 / 2.0) * np.exp(-alpha * beta / p * Rab2)

def kinetic_energy(alpha, RA, beta, RB):
    """Kinetic energy integral for 1s Gaussians."""
    RA = np.array(RA)
    RB = np.array(RB)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    mu = alpha * beta / p
    return mu * (6.0 - 4.0 * mu * Rab2) * (np.pi / p) ** (3.0 / 2.0) * np.exp(-mu * Rab2)

def electron_nucleus_coulomb(alpha, RA, beta, RB, Z, RN):
    """Electron-nucleus Coulomb integral."""
    RA = np.array(RA)
    RB = np.array(RB)
    RN = np.array(RN)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    Rp = (alpha * RA + beta * RB) / p
    Rpn2 = np.dot(Rp - RN, Rp - RN)
    t = p * Rpn2
    return 2 * np.pi * Z / p * np.exp(-alpha * beta / p * Rab2) * F0(t)

def two_electron_coulomb(alpha, RA, beta, RB, gamma, RC, delta, RD):
    """Two-electron Coulomb integral in chemist's notation (ab|cd)."""
    RA = np.array(RA)
    RB = np.array(RB)
    RC = np.array(RC)
    RD = np.array(RD)
    p = alpha + beta
    q = gamma + delta
    pq = p + q
    Rp = (alpha * RA + beta * RB) / p
    Rq = (gamma * RC + delta * RD) / q
    Rpq2 = np.dot(Rp - Rq, Rp - Rq)
    Rac2 = np.dot(RA - RC, RA - RC)
    Rbd2 = np.dot(RB - RD, RB - RD)
    exp_arg1 = -alpha * gamma / p * Rac2
    exp_arg2 = -beta * delta / q * Rbd2
    t = p * q / pq * Rpq2
    return 2 * np.pi ** (5.0 / 2.0) / (p * q * np.sqrt(pq)) * np.exp(exp_arg1 + exp_arg2) * F0(t)

def h2_fci(R, alpha=0.5):
    """
    Solve H2 using Full Configuration Interaction.
    
    Parameters:
    -----------
    R : float
        Bond length (in Bohr)
    alpha : float
        Gaussian exponent for 1s orbitals (default: 0.5)
    
    Returns:
    --------
    energies : array
        Eigenvalues (energies) of the system
    wavefunctions : array
        Eigenvectors (wavefunction coefficients)
    """
    # Nuclear positions
    RA = np.array([0.0, 0.0, -R/2])
    RB = np.array([0.0, 0.0, R/2])
    Z = 1.0  # Nuclear charge
    
    # Compute one-electron integrals
    S_aa = overlap(alpha, RA, alpha, RA)  # = 1.0 for normalized orbitals
    S_ab = overlap(alpha, RA, alpha, RB)
    S_bb = overlap(alpha, RB, alpha, RB)  # = 1.0
    
    T_aa = kinetic_energy(alpha, RA, alpha, RA)
    T_ab = kinetic_energy(alpha, RA, alpha, RB)
    T_bb = kinetic_energy(alpha, RB, alpha, RB)
    
    V_aa_A = electron_nucleus_coulomb(alpha, RA, alpha, RA, Z, RA)
    V_aa_B = electron_nucleus_coulomb(alpha, RA, alpha, RA, Z, RB)
    V_ab_A = electron_nucleus_coulomb(alpha, RA, alpha, RB, Z, RA)
    V_ab_B = electron_nucleus_coulomb(alpha, RA, alpha, RB, Z, RB)
    V_bb_A = electron_nucleus_coulomb(alpha, RB, alpha, RB, Z, RA)
    V_bb_B = electron_nucleus_coulomb(alpha, RB, alpha, RB, Z, RB)
    
    # One-electron Hamiltonian matrix elements
    h_aa = T_aa + V_aa_A + V_aa_B
    h_ab = T_ab + V_ab_A + V_ab_B
    h_bb = T_bb + V_bb_A + V_bb_B
    
    # Two-electron integrals (chemist's notation)
    # (aa|aa), (aa|ab), (aa|bb), (ab|ab), (ab|bb), (bb|bb)
    I_aaaa = two_electron_coulomb(alpha, RA, alpha, RA, alpha, RA, alpha, RA)
    I_aaab = two_electron_coulomb(alpha, RA, alpha, RA, alpha, RA, alpha, RB)
    I_aabb = two_electron_coulomb(alpha, RA, alpha, RA, alpha, RB, alpha, RB)
    I_abab = two_electron_coulomb(alpha, RA, alpha, RB, alpha, RA, alpha, RB)
    I_abbb = two_electron_coulomb(alpha, RA, alpha, RB, alpha, RB, alpha, RB)
    I_bbbb = two_electron_coulomb(alpha, RB, alpha, RB, alpha, RB, alpha, RB)
    
    # Nuclear-nuclear repulsion
    V_nuc = 1.0 / R
    
    # Build 6x6 Hamiltonian matrix
    # Basis: |aa⟩, |bb⟩, |ab↑↑⟩, |ab↓↓⟩, |ab↑↓⟩, |ab↓↑⟩
    H = np.zeros((6, 6))
    
    # Diagonal elements
    # |aa⟩: both electrons on atom A
    H[0, 0] = 2*h_aa + I_aaaa + V_nuc
    
    # |bb⟩: both electrons on atom B
    H[1, 1] = 2*h_bb + I_bbbb + V_nuc
    
    # |ab↑↑⟩: both spin up, one on each atom
    H[2, 2] = h_aa + h_bb + I_aabb + V_nuc
    
    # |ab↓↓⟩: both spin down, one on each atom
    H[3, 3] = h_aa + h_bb + I_aabb + V_nuc
    
    # |ab↑↓⟩: opposite spins, one on each atom
    H[4, 4] = h_aa + h_bb + I_aabb - I_abab + V_nuc
    
    # |ab↓↑⟩: opposite spins, one on each atom (different ordering)
    H[5, 5] = h_aa + h_bb + I_aabb - I_abab + V_nuc
    
    # Off-diagonal elements
    # |aa⟩ <-> |bb⟩: differ by 2 orbitals
    H[0, 1] = I_aaab + I_abbb
    H[1, 0] = H[0, 1]
    
    # |aa⟩ <-> |ab↑↓⟩ and |ab↓↑⟩: differ by 1 orbital
    H[0, 4] = h_ab
    H[4, 0] = h_ab
    H[0, 5] = h_ab
    H[5, 0] = h_ab
    
    # |bb⟩ <-> |ab↑↓⟩ and |ab↓↑⟩: differ by 1 orbital
    H[1, 4] = h_ab
    H[4, 1] = h_ab
    H[1, 5] = h_ab
    H[5, 1] = h_ab
    
    # |ab↑↓⟩ <-> |ab↓↑⟩: differ by 2 orbitals (exchange)
    H[4, 5] = -I_abab
    H[5, 4] = -I_abab
    
    # Solve the eigenvalue problem
    energies, wavefunctions = eigh(H)
    
    return energies, wavefunctions

# Example: Solve H2 at different bond lengths
R_values = [0.5, 1.0, 1.4, 2.0, 3.0]
alpha = 0.5

print("H2 FCI Results (alpha = {:.2f})".format(alpha))
print("=" * 60)
print("{:>8} {:>20} {:>20}".format("R (Bohr)", "E_0 (Hartree)", "E_0 + 1/R"))
print("-" * 60)

for R in R_values:
    energies, _ = h2_fci(R, alpha)
    E0 = energies[0]
    print("{:>8.2f} {:>20.10f} {:>20.10f}".format(R, E0, E0 + 1.0/R))

print("\n")
print("Note: E_0 + 1/R gives the total energy including nuclear repulsion.")
print("The minimum should occur near the equilibrium bond length (~1.4 Bohr).")

Physical Interpretation

The FCI solution for H2 captures:

Covalent bonding: The ground state is dominated by configurations where electrons are shared between atoms
Ionic contributions: Configurations |aa⟩ and |bb⟩ represent ionic states (H⁻H⁺)
Electron correlation: The linear combination of all configurations captures electron-electron correlation effects
Spin coupling: The singlet ground state is a combination of configurations with proper spin symmetry

Comparison with Hartree-Fock

Hartree-Fock (HF) uses only a single Slater determinant, while FCI includes all possible configurations. For H2:

HF: Uses only the bonding configuration, missing ionic and correlation effects
FCI: Includes all 6 configurations, giving the exact solution within the basis set
Energy difference: The difference between HF and FCI energies is the correlation energy

Bond Dissociation

As the bond length $R \to \infty$ , H2 dissociates into two H atoms. The FCI wavefunction correctly describes this:

At large $R$ , the ionic configurations (|aa⟩ and |bb⟩) become less important
The ground state approaches a 50-50 mixture of |ab↑↓⟩ and |ab↓↑⟩
This gives the correct dissociation limit: two separate H atoms

Hartree-Fock fails to describe dissociation correctly because it cannot represent the proper spin coupling at large distances.

References

For more details on FCI, see the Full Configuration Interaction page.

For the integral formulas used here, see the 1s Integrals page.