“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Hartree-Fock Method

Quantum Chemistry

⛓ What you need to know first 6 concepts, 3 layers

The requisite-knowledge inventory for this page, bottom-up: the primitives at the base, combined upward until you reach what this page assumes. Skim the layers you already own; start wherever the ground gets unfamiliar.

base
- Linear algebraconcept
- Quantum mechanics (states & operators)
L1
L2
- Variational principle
↳you are here

1 of these are concepts without a dedicated page yet — the grey chips. Following the linked ones first makes the rest land.

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.

Know this cold

Essential facts

The energy HF misses is the correlation energy, defined as exactly what is left over from the exact result.
$E_{corr} = E_{exact} - E_{HF} < 0$
Only a percent or two of the total energy, but it dominates the energy differences chemistry lives on — bond strengths, barriers.
where this comes from →
Exchange ( $\hat{K}$ ) has no classical analog and acts only between same-spin electrons. It is the energy of the Fermi hole — the avoidance antisymmetry forces on parallel spins.
Opposite-spin electrons feel only Coulomb, no exchange — which is why HF misses opposite-spin correlation entirely.
where this comes from →
The Fock operator is the one-electron core ( $\hat{h}$ ) plus the averaged Coulomb minus exchange of every occupied orbital.
$\hat{F} = \hat{h} + b \sum (\hat{J}_{b} - \hat{K}_{b})$
Its eigenfunctions are the molecular orbitals; its eigenvalues are the orbital energies $ε_{i}$ .
where this comes from →
Koopmans' theorem: minus an occupied orbital energy approximates the ionization potential, assuming the remaining orbitals don't relax.
$IP \approx - ε_{i}$
It works better than it should by error cancellation — neglected orbital relaxation and neglected correlation pull in opposite directions.
where this comes from →
Each electron moves in the average field of all the others, not their instantaneous positions. The coupled $N$ -electron problem becomes $N$ one-electron problems solved self-consistently.
Replacing the instantaneous repulsion with an average is what makes HF solvable — and exactly what it gets wrong (correlation).
where this comes from →
Restricted HF forces both electrons of a pair into the same spatial orbital. This breaks at dissociation: RHF cannot separate H $_{2}$ into two neutral atoms, and the energy sails high as the bond stretches.
The fix (UHF) lets the spins localize on different atoms, but contaminates the spin state — broken symmetry, $⟨ S^{2} ⟩ \neq = 0$ .
where this comes from →
HF is solved by iteration, not in one shot: the Fock operator is built from the very orbitals you are solving for. Guess the density, build the field, solve, repeat until the orbitals stop changing.
$\hat{F} [{ϕ_{i}}] ϕ_{i} = ε_{i} ϕ_{i}$
It is a fixed-point iteration — "self-consistent" means the field reproduces the orbitals that produced it.
where this comes from →
Hartree–Fock approximates the full many-electron wavefunction by a single Slater determinant — one antisymmetrized configuration. Every other choice in the method is just picking the best orbitals to fill it.
This one-determinant ansatz is the source of everything HF gets wrong: the true wavefunction is a sum of many determinants.
where this comes from →

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.

From AO integrals to MO integrals

Every post-HF method listed below — CI, CC, MP, TDDFT — needs the two-electron integrals in the molecular-orbital (MO) basis, not the atomic-orbital (AO) basis you started with. Two different bases for the same operators:

AO basis ${ϕ_{μ}}$ : the Gaussians centered on each atom. This is what you choose at the start of a calculation (6-31G, cc-pVDZ, etc.) and what the integral evaluator actually computes over. Cheap to evaluate, physically arbitrary.
MO basis ${ψ_{i}}$ : linear combinations of AOs whose coefficients $C_{μ i}$ come from diagonalizing the converged Fock matrix. Physically meaningful — occupied MOs hold the electrons, virtual MOs are the unoccupied states.

The two bases are related by the coefficient matrix produced by SCF:

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

Post-HF methods need quantities like "the exchange interaction between MO $i$ and MO $a$ " — indexed by MOs, not AOs. So after SCF converges, you do an integral transformation: rotate the four-index electron-repulsion tensor from AO indices to MO indices.

(pq ∣ rs)_{MO} = μν λσ \sum C_{μ p} C_{ν q} C_{λ r} C_{σ s} (μν ∣ λσ)_{AO}

Done naively as one big sum this is $O (N^{8})$ . Done as four successive single-index contractions — first $μ \to p$ , then $ν \to q$ , then $λ \to r$ , then $σ \to s$ — it's four $O (N^{5})$ steps. Every serious quantum chemistry package implements this transformation; in PySCF it's a single call: ao2mo.incore.full(eri_ao, mf.mo_coeff).

The same idea applies to one-electron operators like $\hat{h}$ : the AO-basis $H_{μν}^{core}$ becomes the MO-basis $h_{ij} = \sum_{μν} C_{μ i} C_{ν j} H_{μν}^{core}$ — same matrix, expressed in a different basis. The Total Energy section above already used this for the diagonal entries $h_{ii}$ .

Once the MO-basis integrals are in hand, every post-HF method below consumes them as input. The Casida equation for TDHF/TDDFT, for example, builds its A and B matrices from a specific subset of these (occupied-virtual orbital pairs) plus the orbital energies $ε_{i}$ from the SCF diagonalization.

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

Q&A

Study notes for the page above. Framing first, then drill-down, then a final compression. Pair them with the derivation when you want to test recall.

What problem does Hartree-Fock solve, and why is the many-body Schrödinger equation otherwise intractable?

The electronic structure of any molecule. $\hat{H} Ψ = E Ψ$ for $N$ electrons is a $3 N$ -dimensional PDE — for water (10 electrons), that's 30 spatial coordinates. A direct numerical solution on a grid of 100 points per dimension would be $10 0^{30} = 1 0^{60}$ grid points, which is hopeless. Hartree-Fock approximates the $N$ -electron wavefunction as a single Slater determinant of $N$ one-electron orbitals, collapsing the problem to $N$ coupled 3D equations — one per orbital. The cost is that HF ignores any correlation between electrons beyond what antisymmetry already enforces.

What's the core ansatz?

The $N$ -electron wavefunction is a single Slater determinant of $N$ spin-orbitals:

Ψ_{HF} = \frac{1}{N !} det [χ_{1} (1) χ_{2} (2) \dots χ_{N} (N)] .

Each $χ_{i}$ is a single-particle wavefunction (orbital times spin). The determinant structure forces $Ψ$ to be antisymmetric under particle exchange, satisfying the Pauli principle. The free parameters are the orbitals; you find them by minimizing $⟨ Ψ_{HF} ∣ \hat{H} ∣ Ψ_{HF} ⟩$ .

Why a determinant and not just a product of orbitals?

A simple Hartree product $Ψ (1, 2, \dots, N) = ψ_{1} (1) ψ_{2} (2) \dots$ is not antisymmetric — swap particles 1 and 2 and you get a different function, not the negative of the original. Fermions require $\hat{P}_{ij} Ψ = - Ψ$ . The determinant builds that in by linearly combining all $N!$ permutations with alternating signs. Hartree's original method used the product without antisymmetrization — it gives wrong results because it misses exchange.

How do the Fock equations come out of the variational principle?

Minimize $⟨ Ψ_{HF} ∣ \hat{H} ∣ Ψ_{HF} ⟩$ over the orbitals $χ_{i}$ subject to the orthonormality constraint $⟨ χ_{i} ∣ χ_{j} ⟩ = δ_{ij}$ . Use Lagrange multipliers (one per orthonormality constraint). The Euler-Lagrange equations are $\hat{F} χ_{i} = ε_{i} χ_{i}$ — a one-particle eigenvalue problem for each orbital. The $ε_{i}$ are the Lagrange multipliers (diagonalized by a canonical-orbital transformation) and we read them as orbital energies.

What's the Fock operator made of?

Three pieces:

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

$\hat{h}$ is the one-electron operator (kinetic energy plus nuclear attraction) — the same as in the bare Schrödinger equation. $\hat{J}_{j}$ is the Coulomb operator for orbital $j$ : the average Coulomb potential from orbital $j$ 's charge density. $\hat{K}_{j}$ is the exchange operator: a nonlocal contribution that captures antisymmetry-induced same-spin avoidance. The $2 J - K$ structure reflects spin: each spatial orbital is doubly occupied (two electrons, one of each spin), so the Coulomb sum is doubled, but exchange only acts between same-spin pairs.

What are J and K physically?

$\hat{J}_{j}$ is classical: the Coulomb potential from the charge density of orbital $j$ , $\hat{J}_{j} (1) = \int ∣ χ_{j} (2) ∣^{2} / r_{12} d 2$ . Local in position — acts as multiplication by a potential $V_{j} (r_{1})$ . $\hat{K}_{j}$ is quantum-mechanical: the exchange-hole effect from antisymmetry, $\hat{K}_{j} (1) χ_{i} (1) = [\int χ_{j}^{*} (2) χ_{i} (2) / r_{12} d 2] χ_{j} (1)$ . Nonlocal — same-spin electrons can't be at the same point, so they effectively avoid each other beyond Coulomb. $\hat{K}$ lives in HF and hybrid DFT; pure DFT functionals approximate it with a local exchange-correlation potential.

Why is SCF iterative?

The Fock operator $\hat{F}$ depends on the orbitals $χ_{i}$ (through $\hat{J}$ and $\hat{K}$ , which sum over occupied orbitals). So $\hat{F} χ_{i} = ε_{i} χ_{i}$ isn't a standard eigenvalue problem — the operator changes when its eigenvectors do. Fix: start with a guess for the orbitals (often from an extended Hückel or a previous calculation), build $\hat{F}$ from them, diagonalize to get new orbitals, build $\hat{F}$ again, repeat until self-consistent. Typically 10-30 iterations for a closed-shell molecule near equilibrium; failures are common in difficult cases (open shells, near-degeneracy, broken symmetry).

What's the difference between RHF and UHF?

RHF (Restricted HF) constrains each spatial orbital to be doubly occupied — one spin-up and one spin-down electron share the same spatial wavefunction. Half the variables, half the orbitals to diagonalize. Works for closed-shell systems near equilibrium. UHF (Unrestricted HF) lets the spin-up and spin-down spatial orbitals differ — $α$ and $β$ orbitals are computed independently. Twice the variables, but can describe open-shell systems and broken-symmetry states. The cost: UHF wavefunctions can be spin-contaminated (a "singlet" UHF result can have triplet character). The right choice depends on what physics you need.

What's Koopmans' theorem?

$- ε_{HOMO} \approx$ ionization potential (IP), and $- ε_{LUMO} \approx$ electron affinity (EA). The approximation: ionizing or attaching an electron doesn't relax the other orbitals, so the energy cost is just removing or adding an electron from a specific orbital — equal to that orbital's energy with a sign flip. Works surprisingly well for closed-shell molecules (errors ~ 1 eV), breaks down for transition metals, open-shell species, and any system where orbital relaxation upon ionization is large.

What's the correlation energy and why does HF miss it?

$E_{correlation} = E_{exact} - E_{HF}$ . HF treats electron-electron repulsion through an AVERAGE field (the $J - K$ combination), where each electron sees the mean potential from all others. This misses the INSTANTANEOUS correlation between electron positions — electrons actually correlate their motion to avoid each other beyond what antisymmetry enforces. For H₂ near equilibrium, correlation energy is ~1% of the total. For stretched bonds, where correlation is "static" (two near-degenerate configurations contribute), HF can be qualitatively wrong by 50% or more.

Why does HF fail at bond dissociation?

At a stretched H₂ bond, the true wavefunction wants to be a roughly equal mixture of two configurations: $σ_{g}^{2}$ (bonding) and $σ_{u}^{2}$ (antibonding). A single Slater determinant can only be ONE of these — RHF stays as $σ_{g}^{2}$ , which puts both electrons on the same atom 50% of the time at infinite separation, giving an ionic dissociation product (H^- + H⁺) instead of two neutral H atoms. Past the Coulson-Fischer point UHF can recover the correct energy by breaking spin symmetry, but at the cost of spin contamination. The proper fix is multireference methods (CASSCF) that mix the two determinants explicitly.

What's the cleanest mental compression?

Three steps. (1) Pick a single Slater determinant as your wavefunction — $N$ spin-orbitals, antisymmetric by construction. (2) Minimize $⟨ Ψ∣ \hat{H} ∣Ψ ⟩$ over the orbitals, subject to orthonormality. (3) The Euler-Lagrange equations are the Fock equations $\hat{F} χ_{i} = ε_{i} χ_{i}$ , where $\hat{F}$ depends on the orbitals through Coulomb and exchange, making the problem self-consistent. Everything else — RHF vs UHF, Koopmans, correlation energy, bond-dissociation failure — is a consequence of the single-determinant ansatz.

References

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

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

Standard references for everything derived on this page: A. Szabo and N. S. Ostlund, Modern Quantum Chemistry (Dover, 1996) — the Slater determinant, Fock operator, Slater–Condon rules, Roothaan equations, and the H₂/STO-3G benchmark all follow its conventions; T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (Wiley, 2000) for the second-quantized formulation; C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951) for the matrix form of the SCF equations.