Hartree-Fock and the self-consistent field

Hartree-Fock does not approximate the Schrödinger equation; it restricts the search. The trial state is a single Slater determinant of N orthonormal spin-orbitals, and the energy is minimised over those orbitals themselves — an infinite-dimensional variation, not a linear one. Because the trial state is normalised and antisymmetric, the result is a rigorous upper bound: E_HF ≥ E_exact, always. Slater–Condon rules reduce the expectation value to integrals over one and two orbitals at a time, E = Σᵢ hᵢᵢ + ½ ΣᵢΣⱼ (Jᵢⱼ − Kᵢⱼ), where hᵢᵢ is the kinetic plus nuclear-attraction integral, Jᵢⱼ = ⟨ij|ij⟩ is the classical repulsion between the charge clouds |χᵢ|² and |χⱼ|², and Kᵢⱼ = ⟨ij|ji⟩ is an integral with no classical analogue that survives only when χᵢ and χⱼ carry the same spin. Everything downstream — the operator, the loop, the missing correlation — follows from having chosen one determinant and nothing more.

Minimise E subject to ⟨χᵢ|χⱼ⟩ = δᵢⱼ. Each constraint buys a Lagrange multiplier εᵢⱼ, and stationarity gives f̂χᵢ = Σⱼ εᵢⱼ χⱼ with f̂ = ĥ + Σⱼ (Ĵⱼ − K̂ⱼ). The multiplier matrix is Hermitian, so some unitary transformation diagonalises it — and that same unitary mixing of the occupied orbitals leaves the determinant unchanged up to a phase, so it moves neither the density nor the energy. In that canonical gauge the condition becomes an ordinary eigenvalue problem, f̂χᵢ = εᵢχᵢ, which is the only reason orbital energies exist at all. Read the two-electron pieces carefully. Ĵⱼψ(r) = [∫ |χⱼ(r′)|²/|r − r′| dr′] ψ(r) multiplies ψ by a potential. K̂ⱼψ(r) = [∫ χⱼ*(r′)ψ(r′)/|r − r′| dr′] χⱼ(r) carries the value of ψ at r′ over to r: a nonlocal integral operator, not a potential at all. Put j = i and the two cancel exactly, so an electron never repels itself. Hartree's earlier scheme, with no K̂, let it.

f̂ is assembled from the orbitals it is meant to produce. That circularity is what makes the equations nonlinear, and there is no direct solution — only a fixed-point loop. Guess a density, build f̂ from it, diagonalise, occupy the N lowest orbitals by aufbau, form a new density, and repeat until the density stops moving; convergence is usually declared when the RMS change in the density matrix falls below about 10⁻⁸ and the energy below 10⁻¹⁰ hartree. The loop is not guaranteed to converge. Linearise about the fixed point and each cycle multiplies the density error by a factor λ, so the error goes as λⁿ and the loop converges only while |λ| < 1; for λ < −1 the error flips sign every cycle and grows, which in atoms shows up as charge sloshing between shells, while λ = −1 sits on the boundary and oscillates forever without decaying. Damping fixes most of it: feed forward P_new = (1 − α)P_old + αP_calc, which replaces λ by r = (1 − α) + αλ, so for λ < 1 the damped loop converges whenever 0 < α < 2/(1 − λ). Production codes use DIIS instead, extrapolating from several stored Fock matrices, which routinely turns a hundred cycles into fifteen.

Numerically you do not solve for functions, you expand them. Write χᵢ = Σ_μ C_μi φ_μ over M fixed basis functions — Gaussians centred on nuclei, in practice — and the integro-differential equation collapses to matrices: FC = SCε, with F_μν = ⟨φ_μ|f̂|φ_ν⟩ and S_μν = ⟨φ_μ|φ_ν⟩. S is not the identity, because atom-centred functions overlap, so this is a generalised eigenvalue problem. Hand it straight to scipy.linalg.eigh(F, S), or orthogonalise first with X = S^(−1/2), solve X†FX C′ = C′ε, and recover C = XC′. The cost is lopsided. Building F needs the two-electron integrals (μν|λσ), about M⁴/8 of them once permutational symmetry is used — 7.8 × 10⁵ for a modest M = 50 — while the diagonalisation inside the loop is only O(M³). That is why integrals are computed once and stored, or recomputed on the fly in direct SCF when memory runs short, and why the basis size rather than the cycle count sets the bill.

Three things are easy to misread. First, E ≠ Σᵢ εᵢ. Each εᵢ is hᵢᵢ + Σⱼ (Jᵢⱼ − Kᵢⱼ), so summing over occupied orbitals counts every pair twice, and the total is E = Σᵢ εᵢ − ½ ΣᵢΣⱼ (Jᵢⱼ − Kᵢⱼ). Second, Koopmans' theorem, IEᵢ ≈ −εᵢ, holds only if the other N − 1 orbitals stay frozen and neither state carries correlation. Both assumptions fail, and in opposite directions: relaxation lowers the ion and so lowers the true ionisation energy, while correlation stabilises the neutral more than the ion and raises it. In helium −ε₁ₛ = 24.98 eV against 24.587 eV measured, because +1.53 eV of relaxation error nearly cancels −1.14 eV of correlation error. Third, a virtual orbital is not an excited state: it is an eigenfunction of a Fock operator built from N electrons, so it feels N of them rather than N − 1 and comes out too high and too diffuse. Where accuracy matters, take the difference of two separately converged calculations — ΔSCF — not an eigenvalue.

Antisymmetry alone keeps same-spin electrons apart: exchange digs a Fermi hole around each electron, and Hartree-Fock carries it exactly. What one determinant cannot represent is the Coulomb hole — the instantaneous avoidance of electrons regardless of spin — because each electron is assigned an orbital that does not depend on where the others happen to be. Löwdin's definition names the deficit, E_corr = E_exact − E_HF: −0.042 Ha (−1.14 eV) in helium, −0.094 Ha in beryllium, −0.390 Ha (−10.6 eV) in neon. As a fraction of the total that is 1.4% for helium, 0.6% for beryllium and 0.3% for neon — negligible-looking, yet each is of the order of a chemical bond energy, which is why Hartree-Fock geometries are usually good and its energetics are not. Brillouin's theorem says where to go next: singly excited determinants have zero matrix element with the converged determinant, so any correction must start from double excitations — which is what MP2, CISD and coupled cluster do. One determinant also fails qualitatively wherever a second is nearly degenerate with it: a stretched bond, or beryllium's 2s and 2p.