The Roothaan-Hall equations

Quantum Chemistry

The previous page wrote the Hartree-Fock eigenvalue problem as $\hat{F} χ_{i} = ε_{i} χ_{i}$ , but that's an equation in continuous-function space — infinitely many degrees of freedom for the orbital, infinitely many possible eigenstates of $\hat{F}$ . We need to project it onto a finite basis to put it on a computer. Doing so turns it into a matrix equation. The Roothaan-Hall equations are that matrix form, and they have a small twist that's worth pulling apart: the basis is non-orthogonal, so what would have been a plain eigenvalue problem becomes a generalized one.

From operator to matrix

Expand each unknown spatial orbital in the basis ${g_{i}}$ :

ϕ_{a} (r) = i \sum C_{ia} g_{i} (r) .

The columns of the coefficient matrix $C \in R^{N \times N}$ are the orbitals; row $i$ tells you how much of basis function $g_{i}$ each orbital contains. Substitute into $\hat{F} ϕ_{a} = ε_{a} ϕ_{a}$ , then project onto each $g_{j}$ by left-multiplying with $⟨ g_{j} ∣ \cdot ⟩$ :

i \sum C_{ia} ⟨ g_{j} ∣ \hat{F} ∣ g_{i} ⟩ = ε_{a} i \sum C_{ia} ⟨ g_{j} ∣ g_{i} ⟩ .

Define the Fock matrix $F_{ji} = ⟨ g_{j} ∣ \hat{F} ∣ g_{i} ⟩$ and the overlap matrix $S_{ji} = ⟨ g_{j} ∣ g_{i} ⟩$ . The projection is the matrix equation

F C = S C ε

where $ε$ is the diagonal matrix of orbital energies. This is a generalized eigenvalue problem: a regular eigenvalue problem $A x = λ x$ on the right-hand side has the identity, here it has $S$ . The reason is that our basis is not orthonormal — the Gaussian primitives have nonzero overlap with each other.

The orthonormality constraint on orbitals translates to $C^{T} S C = I$ : the columns of $C$ are " $S$ -orthonormal." This is the discrete image of $⟨ ϕ_{a} ∣ ϕ_{b} ⟩ = δ_{ab}$ — a standard fact about non-orthogonal bases that anyone who's used Gram-Schmidt has met before.

One eigenvalue equation, many orbitals

Diagonalizing an $N \times N$ matrix gives $N$ eigenvalues and eigenvectors at once. The Roothaan-Hall solve therefore yields every orbital the basis can express, not just the occupied ones. With $N_{occ}$ electrons in $N_{occ} /2$ doubly occupied spatial orbitals, we keep the $N_{occ} /2$ lowest as occupied; the rest are the virtual orbitals. Virtuals don't host any electrons in HF, but they're not waste — they're the discrete image of the unoccupied spectrum and they're what every post-HF correlation method starts from.

For helium with three primitive Gaussians, the basis dimension is three. Diagonalization gives one occupied orbital and two virtuals. The occupied orbital is the analogue of 1s; the two virtuals are crude images of higher spherically-symmetric excitations.

The code (new this page)

scipy.linalg.eigh(F, S) handles the generalized Hermitian case directly. We wrap it in solve_roothaan for clarity, and demo it on the helium core Hamiltonian — the part of $\hat{F}$ with the electron-electron repulsion turned off. That's the He⁺ problem: a single electron in the field of a $Z = 2$ nucleus, no other electrons present.

from scipy.linalg import eigh


def solve_roothaan(F, S):
    """Solve the generalized eigenvalue problem  F C = S C diag(eps).

    scipy.linalg.eigh handles the Hermitian generalized case directly,
    returning eigenvalues sorted ascending and basis-coefficient
    eigenvectors that satisfy  C^T S C = I.
    """
    eps, C = eigh(F, S)
    return eps, C


def helium_core_orbitals():
    """Solve the bare core Hamiltonian (no electron-electron repulsion).

    This is the He+ problem in the STO-3G primitive basis: a single
    electron moving in the field of a Z=2 nucleus. Eigenvalues will
    be the hydrogenic-ion 1s, 2s, 3s approximated in a tiny 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

    return solve_roothaan(h, S)

Solving with just the core Hamiltonian isn't HF — there's no electron-electron interaction in $h$ alone. But it gives us a useful baseline. The exact ground-state energy of He⁺ is $- Z^{2} /2 = - 2$ Hartree; the lowest-eigenvalue STO-3G estimate should land close to that, with higher eigenvalues representing crude approximations to excited He⁺ states. We'll also need the core Hamiltonian eigenvectors as the starting guess for the SCF iteration on the next page.

Run it

He+ orbitals (core Hamiltonian only, STO-3G primitive basis):
  orbital 0   eps = -1.968656   (compare exact -2.000000)
  orbital 1   eps = -0.127134
  orbital 2   eps = +6.603892

The lowest STO-3G eigenvalue is $- 1.969$ Hartree, about $31$ millihartree above the analytic He⁺ ground-state energy of $- 2$ . That gap is purely basis incompleteness — three $s$ -type primitives can't fully capture the cusp at the nucleus. The second eigenvalue is the basis's best attempt at a higher hydrogenic state and is much less accurate; the third is a high-energy artefact. This is the standard pattern for any finite-basis quantum chemistry calculation: the lowest orbital is described decently, then the quality deteriorates rapidly.

Now we have three of the four ingredients HF needs: integrals over the basis, a Fock-matrix builder, and a way to solve $F C = S C ε$ . The missing piece is the iteration that ties them together — the Fock matrix depends on $C$ , and we need a procedure that finds a self-consistent $C$ . That's the next page.

Why is the basis non-orthogonal in the first place?

Because we chose Gaussian primitives centred on the same nucleus (later: on different nuclei), and same-centred or close-centred Gaussians overlap. Orthogonalising them by hand is possible — Löwdin or Gram-Schmidt — but it just adds work without changing the answer; scipy.linalg.eigh handles the generalized problem directly. The deeper reason for living with non-orthogonality is pragmatic: chemistry-relevant basis functions (atom-centred, with the right cusp shape) are naturally non-orthogonal between nearby atoms, and forcing orthogonalisation on top of them produces functions with unphysical oscillations that hurt convergence with basis size.

What's the relationship between the orbital energies and the total energy?

The orbital energies are eigenvalues of $\hat{F}$ ; the total energy is the expectation value of the full $\hat{H}$ on the Slater determinant. They differ by the double-counting term we worked out on the previous page: $E = \sum_{a} ε_{a}$ over occupied orbitals minus the over-counted Coulomb-minus-half-exchange term. So summing orbital energies is wrong by a known correction; the correction is structurally important and missing it is one of the classic textbook errors.

Why diagonalize both $\hat{F}$ 's occupied subspace AND its virtual subspace at once?

The Roothaan-Hall solve produces all $N$ eigenvectors as a side effect, but in HF only the occupied ones matter for the energy: virtuals don't enter $P$ . The orbital energies of the virtuals individually aren't meaningful (they're sensitive to the basis and the Fock parametrization). What's meaningful is the span of the virtual subspace, because that's the space where post-HF methods build correlated wavefunctions. For a pure HF calculation you could in principle skip the virtuals; we keep them because in three lines of scipy you get them for free, and because the SCF iteration will use the full $C$ matrix.