University Physics V · Multi-Electron Atoms · 12.5
Hartree-Fock & the Self-Consistent Field
Every electron here sees the others only as an averaged cloud, and that cloud is built from the very orbitals you are still solving for. This is how you break the circle — guess, build, diagonalise, repeat — and what the converged answer still leaves out.
Build the model
Connect the measurement to the mechanism.
Hartree-Fock is the variational principle run over the smallest space that already respects antisymmetry: all single Slater determinants. Fix that ansatz and the energy becomes a functional of the N occupied spin-orbitals — E = Σᵢ hᵢᵢ + ½ ΣᵢΣⱼ (Jᵢⱼ − Kᵢⱼ) by the Slater–Condon rules — so minimising it under ⟨χᵢ|χⱼ⟩ = δᵢⱼ is a calculus-of-variations problem with a Hermitian matrix of Lagrange multipliers. Diagonalise that matrix, which costs nothing because a unitary mixing of occupied orbitals leaves the determinant and the density alone, and stationarity becomes an eigenvalue problem: f̂χᵢ = εᵢχᵢ with f̂ = ĥ + Σⱼ (Ĵⱼ − K̂ⱼ).
It looks like a one-electron Schrödinger equation and it is not, because the operator is assembled from the orbitals it is supposed to produce. The equations are nonlinear, so they are solved by driving a guess to a fixed point — the self-consistent field. In a finite basis the whole thing collapses to the Roothaan generalised eigenproblem FC = SCε, one scipy.linalg.eigh call per cycle.
What it costs is correlation: a determinant hands each electron an orbital that does not know where the others are at that instant, so it keeps the Fermi hole exchange gives it free and misses the Coulomb hole entirely — 1.14 eV in helium, 1.4% of that atom's total energy but of the order of a chemical bond in absolute size.
- Simple definition
- Hartree-Fock is the best single Slater determinant: the one whose orbitals minimise the energy, each orbital an eigenfunction of a Fock operator built from the averaged Coulomb and exchange field of all the others.
- Example
- For helium the converged solution gives ε₁ₛ = −0.91796 Ha and EHF = −2.86168 Ha = −77.87 eV — a rigorous upper bound that still sits 1.14 eV above the exact −79.01 eV.
The j = i term is Jᵢᵢ − Kᵢᵢ = 0, so exchange cancels self-interaction exactly — the Hartree method's worst flaw, gone.
hᵢᵢ, Jᵢⱼ, Kᵢⱼ in hartree; i and j run over all N occupied spin-orbitals
One effective one-electron operator carrying the averaged field of all N electrons; only K̂ knows which of them share its spin.
ĥ = −(ħ²/2m)∇² − Ze²/4πε₀r; Ĵⱼ multiplies by a potential, K̂ⱼ is a nonlocal integral operator
Orthonormality enters as a Hermitian multiplier matrix; the unitary rotation that diagonalises it changes nothing physical.
εᵢ in hartree; the boundary condition selects bound orbitals in L²(R³) ⊗ C²
scipy.linalg.eigh(F, S) solves it in one line — but C rebuilds F, so that line sits inside a loop.
M × M matrices over a non-orthogonal basis of φμ; S ≠ I because atom-centred functions overlap
Helium: 2ε₁ₛ = −49.96 eV against EHF = −77.87 eV, a 27.91 eV gap that is one pair repulsion counted twice.
Each εᵢ already carries the repulsion with every other occupied orbital, in hartree
Helium gives −ε₁ₛ = 24.98 eV against 24.587 eV measured, because +1.53 eV of relaxation nearly cancels −1.14 eV of correlation.
Frozen orbitals only: the ion may not relax, and neither state carries correlation
One determinant is the whole ansatz
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: EHF ≥ Eexact, 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.
Varying under orthonormality gives the Fock operator
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.
The operator depends on its own eigenvectors, so you iterate
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 Pnew = (1 − α)Pold + αPcalc, 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.
In a finite basis it becomes FC = SCε
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.
Reading the output: εᵢ, Koopmans, and ΔSCF
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.
The Fermi hole is free; the Coulomb hole is not
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, Ecorr = Eexact − EHF: −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.
Change one variable at a time
Make the relationship visible.
Start at a = 0.50 with k = −1.5: r = −0.25, the deviation flips sign every cycle and takes 14 of them to fall below 1e-8. Push a to 1.00 - plain undamped iteration - and r = −1.5 throws the trace past both guides onto the clip, the charge sloshing that stalls a real SCF. Then set a = 0.40, where r collapses to zero.
AMPLIFICATION r-0.250
MAGNITUDE |r|0.250
DEVIATION SIZE AT n0.0002 x d0
CYCLES TO 1e-8, 999 = capped14 cycles
Live interpretationAMPLIFICATION r: −0.250. MAGNITUDE |r|: 0.250. DEVIATION SIZE AT n: 0.0002 x d0. CYCLES TO 1e-8, 999 = capped: 14 cycles
Catch the common trap
Explain before calculating.
A converged restricted Hartree-Fock calculation on helium returns the orbital energy ε₁ₛ = −0.91796 Ha and a total energy of −2.86168 Ha. Why is the total energy not simply 2ε₁ₛ = −1.83592 Ha?
Choose an answer to test the model.
Practice & worked examples
Reason from the model, then test the result.
EasyA converged restricted Hartree-Fock calculation on helium returns the orbital energy ε₁ₛ = −0.91796 Ha and the Coulomb integral J₁ₛ₁ₛ = 1.02576 Ha for its two occupied spin-orbitals. Find the one-electron integral h₁₁, then the total Hartree-Fock energy, in hartree and in eV. Take 1 Ha = 27.2114 eV.
- Write the orbital energy out in full: ε₁ = h₁₁ + Σⱼ (J₁j − K₁j) over both occupied spin-orbitals. The j = 1 term is J₁₁ − K₁₁ = 0, self-interaction cancelling itself. The j = 2 term is J₁₂ − K₁₂, and K₁₂ = 0 because the two spin-orbitals carry opposite spins. So ε₁ₛ = h₁₁ + J₁ₛ₁ₛ.
- Solve for the one-electron part: h₁₁ = ε₁ₛ − J₁ₛ₁ₛ = −0.91796 − 1.02576 = −1.94372 Ha. That is kinetic energy plus attraction to the Z = 2 nucleus, with no repulsion in it.
- Now the determinant energy: E = Σᵢ hᵢᵢ + ½ ΣᵢΣⱼ (Jᵢⱼ − Kᵢⱼ) = 2h₁₁ + J₁ₛ₁ₛ = −3.88744 + 1.02576 = −2.86168 Ha.
- Cross-check by the other route: E = Σᵢ εᵢ − ½ ΣᵢΣⱼ (Jᵢⱼ − Kᵢⱼ) = 2(−0.91796) − 1.02576 = −1.83592 − 1.02576 = −2.86168 Ha. Same number, as it must be.
- Convert: −2.86168 × 27.2114 = −77.87 eV.
Answerh₁₁ = −1.94372 Ha; EHF = −2.86168 Ha = −77.87 eV. Adding the two orbital energies instead would have given −49.96 eV, wrong by 27.91 eV — the pair repulsion counted twice.
MediumAt the Hartree-Fock limit helium has EHF = −2.86168 Ha and ε₁ₛ = −0.91796 Ha. He⁺ has one electron, so Hartree-Fock is exact for it: E = −2.00000 Ha. The exact non-relativistic energy of helium is −2.90372 Ha. Find the Koopmans and ΔSCF ionisation energies, the orbital relaxation energy and the correlation energy, and account for the measured 24.587 eV.
- Koopmans: IE = −ε₁ₛ = 0.91796 Ha = 0.91796 × 27.2114 = 24.98 eV. This freezes the ion's orbital at the neutral atom's 1s.
- ΔSCF: IE = EHF(He⁺) − EHF(He) = −2.00000 + 2.86168 = 0.86168 Ha = 23.45 eV. Now each state is separately optimised.
- The difference is the orbital relaxation energy: 0.91796 − 0.86168 = 0.05628 Ha = 1.53 eV. Letting the ion's orbital contract onto the bare Z = 2 nucleus lowers the ion, and so lowers the ionisation energy.
- Correlation: Ecorr(He) = −2.90372 − (−2.86168) = −0.04204 Ha = −1.14 eV. He⁺ has one electron and therefore exactly zero correlation energy, so only the neutral is stabilised and the true ionisation energy goes up.
- Add it back: 23.45 + 1.14 = 24.59 eV, against the measured 24.587 eV.
- Koopmans' 24.98 eV survives because the two neglected effects pull opposite ways: +1.53 eV from frozen orbitals against −1.14 eV from missing correlation, leaving +0.39 eV, or 1.6%.
AnswerKoopmans 24.98 eV, ΔSCF 23.45 eV, relaxation 1.53 eV, correlation −1.14 eV. ΔSCF plus correlation gives 24.59 eV against 24.587 eV measured; Koopmans is the closer estimate here only by cancellation.
HardOne SCF cycle for a model two-electron closed shell in a two-function basis has S = [[1, 0.400], [0.400, 1]], core Hamiltonian H = [[−1.400, −1.100], [−1.100, −1.400]] Ha, and a Fock matrix built from the current density, F = [[−1.200, −0.900], [−0.900, −1.200]] Ha. Solve FC = SCε for the occupied orbital, build the density matrix, find the electronic energy, and compare it with 2ε₁.
- Both basis functions are equivalent, so the eigenvectors are the symmetric and antisymmetric combinations. For c = (1, 1): Fc = (−2.100, −2.100) and Sc = (1.400, 1.400), so ε₊ = −2.100/1.400 = −1.500 Ha. For c = (1, −1): Fc = (−0.300, 0.300) and Sc = (0.600, −0.600), so ε₋ = −0.500 Ha.
- The occupied orbital is the lower root, ε₁ = −1.500 Ha. Normalise it against S, not against the identity: cᵀSc = 1 + 0.400 + 0.400 + 1 = 2.800, so C₁₁ = C₂₁ = 1/√2.800 = 0.5976.
- Closed-shell density: Pμν = 2 Cμ1 Cν1 = 2 × 0.5976² = 0.7143 in every entry. Check the electron count against the overlap: Σμν Pμν Sμν = 0.7143 × 2.800 = 2.00, two electrons.
- Electronic energy: E = ½ Σμν Pνμ (Hμν + Fμν). Here H + F = [[−2.600, −2.000], [−2.000, −2.600]], so Σ P(H + F) = 0.7143 × (−2.600 − 2.000 − 2.000 − 2.600) = 0.7143 × (−9.200) = −6.5714, and E = −3.2857 Ha.
- Compare with 2ε₁ = −3.000 Ha. The 0.2857 Ha gap, 7.77 eV, is the two-electron energy ½ Σ P(F − H) = ½ × 0.7143 × 0.800 = 0.2857 Ha, counted once in E and twice in the orbital energies.
- This is one cycle, not the answer. F was built from the previous density, and the new P differs from it, so the next cycle returns a different ε — and only when P stops moving is the field self-consistent.
Answerε₁ = −1.500 Ha with C₁₁ = C₂₁ = 0.5976; every entry of P is 0.7143; E = −3.2857 Ha, lying 0.2857 Ha (7.77 eV) below 2ε₁ because the orbital energies double-count that much electron repulsion.