Hartree-Fock Method

Quantum Chemistry

The Hartree-Fock (HF) method is a mean-field approximation for solving the electronic Schrödinger equation. It represents the many-electron wavefunction as a single Slater determinant of molecular orbitals, which are optimized to minimize the total energy.

The obstacles to understanding Hartree-Fock

The Hartree-Fock method is difficult to understand because it relies on so many interconnected concepts.It's not really something that is easily definable in a short sentence. It's an object of immense complication with many different ways to understand it. When someone asks you "what is Hartree-Fock?" the best first answer might be "it's complicated".

You first need to decide what you mean by Hartree-Fock. Do you mean the calculation? Do you mean the derivation? The word means many different things and if you've only seen it presented one way, you may get a warped perception. Not to delve into philosophy too much, but learning Hartree-Fock is kind of like X-ray diffraction. You need to break the concept up into a bunch of pieces and look at it from different angles to get an idea of what is happening. I'm going to go over a few common directions you could take.

Calculation
Derivation

Calculation. If you start with the calculation route you immediately are faced with how difficult it is to calculate the basic quantitites in Hartree-Fock. This method is designed for multi-electron systems and the multi-electron part introduces most of the difficulty. In particular, the 6D integrals that need to be calculated. In the past computers calculated things much slower than they do, so some extremely optimized methods were invented to calculate these integrals faster. Then you start studying how to construct gaussian basis functions and now you are completely lost on the core logic of the scf loop. The details of the calculation can cause you to lose the entire picture.

Calculation, continued. At the macro-level the Hartree-Fock calculation is a fixed-point iteration, however, it is a very detailed fixed point iteration. It involves integrals, matrix elements, basis functions, density matrices, and even a tensor contraction, which I would argue is one of the most interesting parts of the method.

For example, the integrals. Szabo and Ostlund use 1s integrals, the simplest form of integrals. They are limited in applicability, but because they are quick to program they are used. The self-consistent field loop contains a number of steps and each step requires the understanding of many different concepts. You read the loop and then you have no idea what each of the components is or what it means. For example, what is a generalized eigenvalue problem? This leads to my final point about concepts. At the highest level you need to understand the loop, but at the lowest level you need to know how the electron-electron repulsion integral is calculated.

Hartree-Fock should be read as single Slater determinant and many-body to single-particle

If you learn anything from this page, remember this: The Hartree-Fock method uses a single Slater determinant and converts a many body problem to a single particle problem.

A few different concepts, one name

In addition to the fact that the Hartree-Fock wavefunction is a single Slater determinant, there are other concepts involved in the Hartree-Fock solution. First off, there is the Hartree product. The Hartree product states that a many-body wavefunction can be written as a product of independent single-particle wavefunctions. Fock and Slater realized that the wavefunction in this form did not fulfill the antisymmetry requirement (the Pauli exclusion principle) for fermions. One can naturally think of the Slater determinant as the logical next step from a Hartree product. Next we need to write the Hamiltonian. The Hamiltonian in the Hartree-Fock case is reorganized in an interesting way. From a physical perspective there are one-electron contributions and two-electron contributions. Now, the Fock operator rewrites these as a sum of three terms. The one-electron piece is a single term. Then the Coulomb operator and the exchange operator. From an application perspective, at least for me, it makes more sense to think of it as a one-electron set of terms plus a two-electron set of terms. Decomposing further into Coulomb and exchange becomes quite confusing for someone trying to get started. Thinking about the problem in terms of the Fock operator, in my opinion, does not do much to help the reader get started. It is a useful dense notation once you understand what is going on, because you can compactly write the problem. However, the densest, most compact notation is not always the most helpful for first-time understanding.

The Hartree-Fock Potential

Hartree-Fock describes an N-electron system, but it does so by solving N coupled single-particle equations, one per orbital. Each orbital sees an effective single-particle potential built from all the others. The Coulomb and exchange terms are single-particle operators. We start with a many-electron Hamiltonian and end up with a single-electron Hamiltonian with an effective potential. The transformation is a calculus-of-variations argument; I'll write up the variational derivation properly when I'm ready to stand behind it.

Slater-Condon Rules

The Slater-Condon rules are extremely important to understanding many-body theory. The rules are the result of a derivation which has a few different forms. From a high level, you take the expectation value of the electronic Hamiltonian over a Slater determinant. You then use the orthogonality of the spin orbitals to reduce the number of terms you have to work with. There are six separate cases to consider when performing the integrations: three associated with the one-electron operator and three associated with the two-electron operator. A confusing bit of nomenclature that you run into is the one-electron operator acting on the many-electron wavefunction. Likewise the two-electron operator acting on the many-electron wavefunction. When you go about deriving the Hartree-Fock equations, you begin by constructing an expectation value of a Slater determinant. Easy enough. However, after this expectation is calculated, by some magical rules (known as the Slater-Condon rules), the determinant has a completely different form. It looks like someone took the determinant of a 2×2 matrix and then somehow found a way to tack a summation sign on front.

Derivation of the Hartree-Fock Equations (functional variation and canonical orbitals)

Basically you apply the variational principle to get the Hartree-Fock equations written down. You separate the electron repulsion integrals into Coulomb and exchange terms for physical interpretation. The only problem is that the equation is not written in an $x = f (x)$ form. There are still two electrons in the notation. In fact, the fixed-point iteration is in some ways a fixed-function iteration. In order to perform the calculation you must propose an initial $ψ_{0}$ from which the self-consistent field method will begin. As an artifact of the variational procedure there is a transformation matrix left in the Hartree-Fock equations. The equation, for a moment, looks like a generalized eigenvalue problem instead of an ordinary eigenvalue problem. This is where one introduces the canonical orbitals. These orbitals are taken to be orthogonal in the Fock basis. From there you see that you can apply a unitary transformation to the Hartree-Fock equation to put it in its final mean-field self-consistent form. The difficulty here is that there are a lot of "technical details" that get you to the final answer: calculus of variations and a unitary transformation.

Wavefunction

The Hartree-Fock wavefunction is a single Slater determinant:

∣ Ψ_{HF} ⟩ = ∣ χ_{1} χ_{2} \dots χ_{N} ⟩

where $χ_{i}$ are spin orbitals (spatial orbitals multiplied by spin functions) and $N$ is the number of electrons.

For a closed-shell system with $N$ electrons, the wavefunction can be written as:

∣ Ψ_{HF} ⟩ = ∣ ψ_{1} \overset{ˉ}{ψ}_{1} ψ_{2} \overset{ˉ}{ψ}_{2} \dots ψ_{N /2} \overset{ˉ}{ψ}_{N /2} ⟩

where $ψ_{i}$ are spatial molecular orbitals and the bar denotes beta (down) spin.

Molecular Orbitals

Molecular orbitals are expanded in a basis set:

ψ_{i} (r) = μ \sum C_{μ i} ϕ_{μ} (r)

where $ϕ_{μ}$ are atomic orbital basis functions and $C_{μ i}$ are the molecular orbital coefficients.

Fock Equation

The Hartree-Fock equations are derived by applying the variational principle to minimize the energy. This leads to the Fock equation:

\hat{F} ∣ ψ_{i} ⟩ = ϵ_{i} ∣ ψ_{i} ⟩

where $\hat{F}$ is the Fock operator and $ϵ_{i}$ are the orbital energies.

In matrix form (using the basis set expansion), this becomes:

FC = SC ϵ

where:

$F$ is the Fock matrix
$C$ is the matrix of molecular orbital coefficients
$S$ is the overlap matrix
$ϵ$ is the diagonal matrix of orbital energies

Fock Operator

The Fock operator is:

\hat{F} (1) = \hat{h} (1) + j = 1 \sum N /2 [2 \hat{J}_{j} (1) - \hat{K}_{j} (1)]

where:

$\hat{h} (1)$ is the one-electron operator (kinetic energy + electron-nucleus attraction)
$\hat{J}_{j} (1)$ is the Coulomb operator
$\hat{K}_{j} (1)$ is the exchange operator

The Coulomb operator is:

\hat{J}_{j} (1) = \int ψ_{j}^{*} (2) \frac{1}{r _{12}} ψ_{j} (2) d r_{2}

and the exchange operator is:

\hat{K}_{j} (1) ψ_{i} (1) = \int ψ_{j}^{*} (2) \frac{1}{r _{12}} ψ_{i} (2) d r_{2} ψ_{j} (1)

Fock Matrix Elements

The Fock matrix elements in the atomic orbital basis are:

F_{μν} = H_{μν}^{core} + λσ \sum P_{λσ} [(μν ∣ λσ) - \frac{1}{2} (μ λ ∣ ν σ)]

where:

$H_{μν}^{core} = ⟨ ϕ_{μ} ∣ \hat{h} ∣ ϕ_{ν} ⟩$ is the core Hamiltonian
$P_{λσ}$ is the density matrix
$(μν ∣ λσ)$ are two-electron integrals in chemist's notation

The density matrix is:

P_{μν} = 2 i = 1 \sum N /2 C_{μ i} C_{ν i}^{*}

where the factor of 2 accounts for the two electrons (alpha and beta) in each occupied orbital.

Self-Consistent Field (SCF) Procedure

The Hartree-Fock equations must be solved iteratively because the Fock operator depends on the molecular orbitals through the density matrix. The SCF procedure is:

Choose an initial guess for the molecular orbital coefficients $C$
Build the density matrix $P$
Construct the Fock matrix $F$
Solve the generalized eigenvalue problem $FC = SC ϵ$
Check for convergence (density matrix or energy)
If not converged, return to step 2 with new orbitals

Total Energy

The Hartree-Fock total energy is:

E_{HF} = μν \sum P_{μν} H_{μν}^{core} + \frac{1}{2} μν λσ \sum P_{μν} P_{λσ} [(μν ∣ λσ) - \frac{1}{2} (μ λ ∣ ν σ)] + V_{nuc}

where $V_{nuc}$ is the nuclear-nuclear repulsion energy.

Alternatively, using orbital energies:

E_{HF} = i = 1 \sum N /2 (ϵ_{i} + h_{ii}) + V_{nuc}

Here $h_{ii} = ⟨ ψ_{i} ∣ \hat{h} ∣ ψ_{i} ⟩$ is the one-electron core integral evaluated in the molecular-orbital basis (not the atomic-orbital $H_{μν}^{core}$ defined above). The two are related by $h_{ii} = \sum_{μν} C_{μ i}^{*} C_{ν i} H_{μν}^{core}$ .

Koopmans' Theorem

Koopmans' theorem states that the negative of the orbital energy is approximately equal to the ionization potential (for occupied orbitals) or electron affinity (for virtual orbitals):

- ϵ_{i} \approx IP_{i}

This approximation neglects orbital relaxation upon electron removal/addition.

Python Implementation

Here is a simplified Python implementation of the Hartree-Fock method for a minimal basis H2 molecule:

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 unnormalized 1s Gaussian primitives.
    Szabo & Ostlund Appendix A: T_ab = mu*(3 - 2*mu*R^2)*(pi/p)^(3/2)*exp(-mu*R^2).
    """
    RA = np.array(RA)
    RB = np.array(RB)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    mu = alpha * beta / p
    return mu * (3.0 - 2.0 * mu * Rab2) * (np.pi / p) ** 1.5 * np.exp(-mu * Rab2)

def electron_nucleus_coulomb(alpha, RA, beta, RB, Z, RN):
    """Electron-nucleus attraction integral for unnormalized 1s Gaussians.
    Szabo & Ostlund Appendix A: V_ab = -(2*pi*Z/p)*exp(-mu*R^2)*F0(p*|RP-RN|^2).
    The minus sign is the attractive electron-nucleus interaction.
    """
    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 repulsion integral (AB|CD) for unnormalized 1s Gaussians.
    Szabo & Ostlund Appendix A:
        (AB|CD) = (2 pi^(5/2) / (p q sqrt(p+q)))
                  * exp(-alpha*beta/p * |RA-RB|^2 - gamma*delta/q * |RC-RD|^2)
                  * F0(p*q/(p+q) * |RP-RQ|^2)
    Electron 1 lives in primitives A, B; electron 2 lives in primitives C, D.
    """
    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)
    Rab2 = np.dot(RA - RB, RA - RB)
    Rcd2 = np.dot(RC - RD, RC - RD)
    K_AB = np.exp(-alpha * beta / p * Rab2)
    K_CD = np.exp(-gamma * delta / q * Rcd2)
    return 2 * np.pi ** 2.5 / (p * q * np.sqrt(pq)) * K_AB * K_CD * F0(p * q / pq * Rpq2)

def hartree_fock_h2(R, alpha=0.5, max_iter=100, tol=1e-8):
    """
    Solve H2 using the Hartree-Fock method.
    
    Parameters:
    -----------
    R : float
        Bond length (in Bohr)
    alpha : float
        Gaussian exponent for 1s orbitals (default: 0.5)
    max_iter : int
        Maximum number of SCF iterations
    tol : float
        Convergence tolerance for density matrix
    
    Returns:
    --------
    energy : float
        Total Hartree-Fock energy
    C : array
        Molecular orbital coefficients
    epsilon : array
        Orbital energies
    """
    # 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
    
    # Number of basis functions (2 for minimal basis H2)
    n_basis = 2
    
    # Compute overlap matrix S
    S = np.zeros((n_basis, n_basis))
    S[0, 0] = overlap(alpha, RA, alpha, RA)  # = 1.0 for normalized
    S[0, 1] = overlap(alpha, RA, alpha, RB)
    S[1, 0] = S[0, 1]
    S[1, 1] = overlap(alpha, RB, alpha, RB)  # = 1.0
    
    # Compute core Hamiltonian H_core
    H_core = np.zeros((n_basis, n_basis))
    
    # One-electron integrals
    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)
    
    H_core[0, 0] = T_aa + V_aa_A + V_aa_B
    H_core[0, 1] = T_ab + V_ab_A + V_ab_B
    H_core[1, 0] = H_core[0, 1]
    H_core[1, 1] = T_bb + V_bb_A + V_bb_B
    
    # Two-electron integrals (chemist's notation)
    # Store as (aa|aa), (aa|ab), (aa|bb), (ab|ab), (ab|bb), (bb|bb)
    # Index: 0=aa, 1=ab, 2=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)
    
    # Store integrals in a 4-index array for convenience
    # I[i,j,k,l] = (ij|kl)
    I = np.zeros((n_basis, n_basis, n_basis, n_basis))
    I[0, 0, 0, 0] = I_aaaa
    I[0, 0, 0, 1] = I_aaab
    I[0, 0, 1, 0] = I_aaab
    I[0, 0, 1, 1] = I_aabb
    I[0, 1, 0, 0] = I_aaab
    I[0, 1, 0, 1] = I_abab
    I[0, 1, 1, 0] = I_abab
    I[0, 1, 1, 1] = I_abbb
    I[1, 0, 0, 0] = I_aaab
    I[1, 0, 0, 1] = I_abab
    I[1, 0, 1, 0] = I_abab
    I[1, 0, 1, 1] = I_abbb
    I[1, 1, 0, 0] = I_aabb
    I[1, 1, 0, 1] = I_abbb
    I[1, 1, 1, 0] = I_abbb
    I[1, 1, 1, 1] = I_bbbb
    
    # Nuclear-nuclear repulsion
    V_nuc = 1.0 / R
    
    # Initial guess: diagonalize the core Hamiltonian in an orthogonal basis.
    # Orthogonalize using S^(-1/2).
    S_eval, S_evec = eigh(S)
    S_inv_half = S_evec @ np.diag(1.0 / np.sqrt(S_eval)) @ S_evec.T

    H_ortho = S_inv_half.T @ H_core @ S_inv_half
    epsilon, C_ortho = eigh(H_ortho)
    C = S_inv_half @ C_ortho

    # Build the initial density matrix from the core-guess orbitals.
    P = np.zeros((n_basis, n_basis))
    for mu in range(n_basis):
        for nu in range(n_basis):
            P[mu, nu] = 2.0 * C[mu, 0] * C[nu, 0]  # only first orbital occupied

    # SCF iteration. We iterate the density matrix P (not the coefficients C)
    # because eigenvectors carry an arbitrary sign that can flip iteration to
    # iteration; mixing P sidesteps that.
    damping = 0.5
    converged = False
    delta_P = float("inf")
    for iteration in range(max_iter):
        # Build Fock matrix from current P.
        F = H_core.copy()
        for mu in range(n_basis):
            for nu in range(n_basis):
                for lam in range(n_basis):
                    for sig in range(n_basis):
                        # Coulomb term: (mu nu | lam sig)
                        # Exchange term: -1/2 (mu lam | nu sig)
                        F[mu, nu] += P[lam, sig] * (
                            I[mu, nu, lam, sig] - 0.5 * I[mu, lam, nu, sig]
                        )

        # Diagonalize F in the orthogonal basis, transform back.
        F_ortho = S_inv_half.T @ F @ S_inv_half
        epsilon, C_ortho = eigh(F_ortho)
        C = S_inv_half @ C_ortho

        # New density from new orbitals.
        P_new = np.zeros((n_basis, n_basis))
        for mu in range(n_basis):
            for nu in range(n_basis):
                P_new[mu, nu] = 2.0 * C[mu, 0] * C[nu, 0]

        delta_P = np.max(np.abs(P_new - P))
        if delta_P < tol:
            P = P_new
            converged = True
            break

        # Density damping for stability.
        P = damping * P + (1.0 - damping) * P_new

    if not converged:
        print(f"Warning: SCF did not converge in {max_iter} iterations "
              f"(final delta_P = {delta_P:.3e})")

    # Compute total energy from the converged density.
    energy = 0.0
    for mu in range(n_basis):
        for nu in range(n_basis):
            energy += P[mu, nu] * H_core[mu, nu]
            for lam in range(n_basis):
                for sig in range(n_basis):
                    energy += 0.5 * P[mu, nu] * P[lam, sig] * (
                        I[mu, nu, lam, sig] - 0.5 * I[mu, lam, nu, sig]
                    )
    energy += V_nuc

    return energy, C, epsilon

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

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

for R in R_values:
    energy, C, epsilon = hartree_fock_h2(R, alpha)
    print("{:>8.2f} {:>20.10f} {:>20.10f}".format(R, energy, energy + 1.0/R))

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

Sanity tests

Integral formulas like the ones above are easy to type wrong, and the SCF loop will happily eat a wrong Hamiltonian and "converge" to numbers that look plausible but aren't physical. The tests below cross-check each integral against an independent calculation: for the one-electron integrals, direct numerical integration in spherical coordinates over a single primitive at the origin; for the two-electron integral, an analytic special case where all four primitives coincide; for the SCF as a whole, a couple of physically-required properties (energy minimum near $R \approx 1.4$ Bohr, HOMO orbital energy in the right ballpark for Koopmans' theorem). They run in milliseconds. Run them whenever you touch the integral code.

import numpy as np

def _numerical_kinetic_self(alpha, R_max=8.0, N=200_000):
    """<exp(-alpha r^2) | -1/2 nabla^2 | exp(-alpha r^2)> by direct integration.
    Spherical symmetry: 4*pi * integral_0^inf r^2 (3*alpha - 2*alpha^2 r^2) exp(-2*alpha r^2) dr.
    """
    rs = np.linspace(1e-6, R_max, N)
    integrand = 4 * np.pi * rs**2 * (3 * alpha - 2 * alpha**2 * rs**2) * np.exp(-2 * alpha * rs**2)
    return np.trapz(integrand, rs)

def _numerical_eN_self(alpha, Z=1.0, R_max=8.0, N=200_000):
    """<exp(-alpha r^2) | -Z/r | exp(-alpha r^2)> = -4*pi*Z * integral_0^inf r exp(-2*alpha r^2) dr."""
    rs = np.linspace(1e-6, R_max, N)
    return np.trapz(-4 * np.pi * Z * rs * np.exp(-2 * alpha * rs**2), rs)

def sanity_tests(verbose=True, tol=1e-3):
    alpha = 0.5
    origin = [0.0, 0.0, 0.0]
    rows = []

    # 1. Kinetic energy of a single unnormalized 1s Gaussian at the origin.
    f = kinetic_energy(alpha, origin, alpha, origin)
    r = _numerical_kinetic_self(alpha)
    rows.append(("kinetic <phi|T|phi>", f, r, abs(f - r) < tol))

    # 2. Electron-nucleus attraction with electron and nucleus at the origin.
    f = electron_nucleus_coulomb(alpha, origin, alpha, origin, 1.0, origin)
    r = _numerical_eN_self(alpha, Z=1.0)
    rows.append(("electron-nucleus <phi|-Z/r|phi>", f, r, abs(f - r) < tol))

    # 3. Two-electron self-integral (00|00) has a closed form for coincident primitives:
    #    (00|00) = pi^(5/2) / (4 * alpha^(5/2))   when all four exponents equal alpha
    #              and all four centers coincide.
    f = two_electron_coulomb(alpha, origin, alpha, origin, alpha, origin, alpha, origin)
    r = np.pi**2.5 / (4 * alpha**2.5)
    rows.append(("two-electron (00|00)", f, r, abs(f - r) < tol))

    # 4. H2 SCF: energy should minimize near R = 1.4 Bohr,
    #    and -eps[HOMO] should be in the few-tenths-of-a-Hartree range
    #    (Koopmans' theorem says it approximates the ionization potential).
    Rs = [0.7, 1.0, 1.4, 2.0, 3.0]
    energies = [hartree_fock_h2(R, alpha)[0] for R in Rs]
    R_min = Rs[int(np.argmin(energies))]
    e_min, _, eps_min = hartree_fock_h2(R_min, alpha)
    rows.append((f"H2 minimum near R~1.4 (got R={R_min})", None, None, abs(R_min - 1.4) <= 0.6))
    rows.append((f"HOMO eps={eps_min[0]:+.3f} bound and reasonable",
                 None, None, -1.0 < eps_min[0] < 0.0))

    if verbose:
        print(f"{'test':<44} {'formula':>14} {'reference':>14}   ok?")
        print("-" * 82)
        for name, f, r, ok in rows:
            f_s = f"{f:>14.6f}" if f is not None else " " * 14
            r_s = f"{r:>14.6f}" if r is not None else " " * 14
            print(f"{name:<44} {f_s} {r_s}   {'PASS' if ok else 'FAIL'}")

    return all(row[3] for row in rows)

# Run them.
print("\n--- sanity tests ---")
ok = sanity_tests()
print("All checks passed." if ok else "FAILED -- investigate before trusting results.")

Expected output: all six rows show PASS. If any row fails, the integral or SCF logic above has been edited into a wrong state — fix that before reading anything else off the H₂ output table. As a candid aside: an earlier version of this page failed three of these (kinetic, electron-nucleus, two-electron) and the SCF still produced plausible numbers anyway. That's exactly the kind of silent failure these tests exist to catch.

Limitations

The Hartree-Fock method has several important limitations:

Electron correlation: HF neglects electron correlation, which is the difference between the exact energy and the HF energy
Bond dissociation: HF fails to describe bond dissociation correctly, especially for systems with multiple bonds
Static correlation: HF cannot describe systems with near-degeneracy effects
Basis set dependence: Results depend on the quality of the basis set

Comparison with FCI

For H2, we can compare Hartree-Fock with Full Configuration Interaction:

HF: Uses a single Slater determinant, missing correlation effects
FCI: Includes all possible configurations, giving the exact solution within the basis set
Correlation energy: $E_{corr} = E_{FCI} - E_{HF}$

For H2 at equilibrium, the correlation energy is typically a few percent of the total energy, but becomes more important at stretched bond lengths.

Extensions

Several methods extend Hartree-Fock to include correlation:

Configuration Interaction (CI): Adds excited configurations to the HF wavefunction
Coupled Cluster (CC): Uses exponential parameterization of the wavefunction
Møller-Plesset Perturbation Theory (MP): Treats correlation as a perturbation
Density Functional Theory (DFT): Uses functionals of the electron density

References

For comparison with FCI, see the H2 Atom and FCI page.

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