The SCF iteration

Quantum Chemistry

The Fock matrix depends on the orbitals, the orbitals are determined by the Fock matrix. Two pages of machinery and we're chasing our tail. The standard way out is to iterate: start with a guess, build $F$ , diagonalize to get new orbitals, build a new $F$ , repeat. If the iteration converges, the fixed-point answer satisfies both equations at once. That fixed point is what "self-consistent field" names, and the iteration that finds it is the SCF cycle.

The fixed-point structure

Write the HF problem as a function from density matrices to density matrices:

P F [P] F eigh (F, S) C occ P^{'} .

A self-consistent solution is a density matrix that maps to itself: $P = F (P)$ . Plain iteration — $P_{k + 1} = F (P_{k})$ — works for well-behaved problems. For atoms it usually converges in a handful of steps; for tougher cases (transition metals, near-degeneracies) the iteration can stall or oscillate, and the standard remedy is DIIS (direct inversion in the iterative subspace), which we won't need.

The starting guess matters less than you might think. Anything reasonable works. The standard choice — and what we'll use — is the eigenvectors of the core Hamiltonian: solve $h C = S C ε$ with no electron-electron repulsion. The lowest occupied orbital from that solve is the He⁺-like orbital; using it to build the first $F$ ignores the missing repulsion, but the iteration adds it back in a few cycles.

The closed-shell iteration

For a closed-shell system with $N_{occ}$ doubly-occupied spatial orbitals, one cycle is:

build $P_{k l} = 2 \sum_{a \in occ} C_{ka} C_{l a}$
build $F = h + J [P] - \frac{1}{2} K [P]$
compute $E = \frac{1}{2} tr [P (h + F)]$
solve $F C = S C ε$
check $∣ E_{k} - E_{k - 1} ∣$ against tolerance; if not converged, go back to (1)

The loop converges when the energy stops moving. Equivalently, one can monitor the change in the density matrix $∥ P_{k} - P_{k - 1} ∥$ — that's a stricter test, but for a small atomic problem either works.

The code (new this page)

A density-matrix helper, the SCF loop itself, and an entry point that runs SCF on helium with the three primitive Gaussians as the AO basis.

def density_matrix(C, n_occ):
    """Closed-shell density matrix from the n_occ lowest-energy orbitals.

        P_kl = 2 * sum_{a in occ} C[k,a] * C[l,a]
    """
    Co = C[:, :n_occ]
    return 2.0 * Co @ Co.T


def scf(h, S, ERI, n_occ, max_iter=50, tol=1e-10, verbose=True):
    """Restricted closed-shell SCF in a fixed AO basis.

    Inputs are the basis-projected one-electron Hamiltonian h, the AO
    overlap S, the four-index ERI tensor, and the number of doubly
    occupied orbitals n_occ. Returns (E_total, eps, C, P) at convergence.
    """
    # initial guess: diagonalize the bare core Hamiltonian
    eps, C = solve_roothaan(h, S)
    P = density_matrix(C, n_occ)

    E_prev = 0.0
    for it in range(max_iter):
        F = build_fock(h, ERI, P)
        E = hf_energy_from_density(P, h, F)

        if verbose:
            print(f"  iter {it:2d}    E = {E:+.10f}    dE = {E - E_prev:+.2e}")

        if abs(E - E_prev) < tol and it > 0:
            return E, eps, C, P

        eps, C = solve_roothaan(F, S)
        P = density_matrix(C, n_occ)
        E_prev = E

    raise RuntimeError("SCF did not converge")


def helium_scf():
    """SCF on the helium atom in the bare 3-primitive STO-3G basis."""
    a = STO_3G_HE_ALPHAS
    n = len(a)
    Z = 2

    S = np.array([[overlap_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    T = np.array([[kinetic_g(a[i], a[j])               for j in range(n)] for i in range(n)])
    V = np.array([[nuclear_attraction_g(a[i], a[j], Z) for j in range(n)] for i in range(n)])
    h   = T + V
    ERI = build_eri_tensor(a)

    return scf(h, S, ERI, n_occ=1)

Notice the change in basis from page 4. There we used the contracted STO-3G orbital — three primitives wired together with hand-supplied coefficients — as the only basis function. Here the basis is the three primitives themselves, with SCF free to choose any linear combination. This gives the iteration room to find a slightly better orbital than the rigid STO-3G fit. We'll see exactly how much better in the output.

Run it

Helium SCF (3-primitive basis, 1 doubly-occupied orbital):
  iter  0    E = -2.7115784567    dE = -2.71e+00
  iter  1    E = -2.8151312634    dE = -1.04e-01
  iter  2    E = -2.8162312450    dE = -1.10e-03
  iter  3    E = -2.8162460833    dE = -1.48e-05
  iter  4    E = -2.8162463049    dE = -2.22e-07
  iter  5    E = -2.8162463082    dE = -3.39e-09
  iter  6    E = -2.8162463083    dE = -5.20e-11

Converged HF energies:
  helium_scf (variational)  : E = -2.816246
  helium_via_fock (fixed d) : E = -2.807784
  difference                : 8.5e-03

Two things to read off this trace. (1) The first iteration starts at $- 2.712$ — too high — because the initial density came from the He⁺ core eigenvector, which describes one electron in a bare $Z = 2$ field. Building $F$ with that density introduces the electron-electron repulsion for the first time, and the iteration responds by relaxing the orbital outward (the second electron pushes against the first). By iteration 3 the energy is within a hundredth of a millihartree of the answer; by iteration 6 it's at the floating-point floor. Convergence is geometric with a small contraction factor — typical for an easy SCF problem.

(2) The converged SCF energy is $- 2.816246$ Hartree, versus $- 2.807784$ for the rigid STO-3G contraction we used on page 4. SCF in the bare-primitive basis is $8.5$ millihartree lower. The STO-3G coefficients were fitted to a Slater orbital, not optimized on this Hamiltonian, so SCF — given access to the same three primitives but freedom to choose any linear combination — finds a meaningfully better orbital. The lesson generalizes: fixed-coefficient calculations on contracted basis sets are upper bounds even within their own primitive subspace; running SCF on the underlying primitives (or in a richer contracted basis) will lower the energy.

What we now have

A working closed-shell HF program. Give it $(h, S, ERI, N_{occ})$ for any system and it returns the variationally optimal HF energy and orbitals. The SCF code is generic — it doesn't care whether the integrals come from one centre or many — so the only thing standing between us and H₂ is the integrals themselves. Pages 8 and 9 build the multi-centre overlap, kinetic, nuclear-attraction, and two-electron integrals that make a molecule possible. Then page 10 runs scf on H₂ and scans the bond length.

Why does the iteration converge at all?

Because for a well-conditioned closed-shell problem, the map $F$ is locally a contraction near the fixed point — small perturbations in $P$ shrink as they pass through one cycle. This isn't a theorem in general (you can construct cases where it fails — open-shell metals, strong correlation), and that's why production codes ship with DIIS, level-shifting, damping, and a small zoo of stabilization tricks. For helium the bare iteration is fine.

Could the SCF transient ever go below the true ground-state energy?

In principle, yes, and that's an important caveat about the variational principle. The bound $E \geq E_{0}$ applies to self-consistent energies — quantities of the form $⟨ ψ ∣ \hat{H} ∣ ψ ⟩$ with $ψ$ normalized. The energy reported during an SCF transient is computed from a density $P$ that came from a previous iteration's Fock, not from the current iteration's. That mismatch can in pathological cases produce values below $E_{0}$ . For helium it doesn't, but you should not rely on the bound mid-iteration; only the converged answer is variationally trustworthy.

Why is the converged energy strictly above the SCF on a richer basis would give?

Because adding more basis functions can only enlarge the variational trial space, and the variational principle then bounds the new minimum from below by the old one. STO-3G has three primitives per atom, all $s$ -type; STO-6G has six; cc-pVDZ adds polarization $p$ -functions; cc-pVTZ adds $p$ and $d$ . Each enlargement contains the previous as a subspace, so each gives a tighter bound. The basis-set limit is the lowest energy a single Slater determinant can produce; everything below it (the correlation energy) is beyond HF's reach no matter how big the basis gets.

Exercises

A full exercise set is available for this topic, structured as one worked example + 7 practice problems (across 7 surface contexts) + 2 pattern-resistant check problems.

Open the Hartree-Fock SCF exercise set → Recognize a problem as a fixed-point iteration — an operator that depends on its own eigenvector — and solve by iterative refinement: linearize with a guess, solve, use the result to rebuild the operator, iterate. Diagnose convergence, oscillation, and the use of damping or DIIS-style acceleration.