The Roothaan-Hall Equations

Every mainstream quantum chemistry package — Gaussian, ORCA, Psi4, Q-Chem, NWChem, PySCF — performs a Roothaan–Hall (or its unrestricted/restricted-open-shell cousin) iteration as the innermost engine of essentially every calculation, whether the final method is Hartree–Fock itself, a post-HF correlated method, or Kohn–Sham density functional theory (whose matrix equations have the identical $FC=SC\varepsilon$ structure with an exchange-correlation potential folded into $F$). Convergence acceleration techniques such as DIIS (direct inversion in the iterative subspace) exist specifically to speed up the Roothaan–Hall fixed-point iteration, since naive alternation between building $F$ from $P$ and diagonalizing can oscillate or converge painfully slowly for large or strongly correlated systems. In drug discovery, Roothaan–Hall-based Hartree–Fock calculations (and their DFT descendants sharing the same matrix machinery) are used to compute partial charges, conformer energies, and reactive intermediates for docking and force-field parameterization pipelines. In materials science, plane-wave and Gaussian-basis periodic codes solve a crystal-momentum-resolved generalization of the same equations to obtain band structures and cohesive energies. In spectroscopy, the same Fock and density matrices feed directly into computing vibrational frequencies (via the Hessian of the electronic energy), NMR shielding tensors, and UV-Vis transition energies through linear-response extensions built on the converged Roothaan–Hall solution. Even machine-learned interatomic potentials, now common in high-throughput materials screening, are frequently trained on energies and forces generated from Roothaan–Hall-type SCF calculations, making this 1951 reformulation one of the most continuously used pieces of infrastructure in all of computational science.

:::quiz[In the closed-shell Fock matrix formula F_μν = H^core_μν + Σ P_λσ[(μν|λσ) − ½(μλ|νσ)], what does the ½ factor multiplying the second term represent?]

• A normalization constant for the basis functions
• The exchange contribution's correct relative weight compared to twice the Coulomb term ✓
• A relativistic correction factor
• An artifact of using atomic units

Ask the tutor: Can you walk through how DIIS actually accelerates convergence of the Roothaan–Hall SCF loop?

Hartree–Fock theory says that each electron in a closed-shell molecule moves in an effective, averaged potential created by the nuclei and by all the other electrons, with exchange included exactly for same-spin pairs. The orbital $\psi_i(\mathbf r)$ that describes this motion satisfies the Fock equation $\hat f\,\psi_i = \varepsilon_i\,\psi_i$, where $\hat f$ is the one-electron Fock operator. This equation looks like an ordinary eigenvalue problem, but it is not one in practice: $\hat f$ contains Coulomb and exchange operators built from the very orbitals we are trying to find, and $\psi_i$ is an arbitrary function of three continuous coordinates, not a vector in a finite space. Roothaan and Hall's insight was to stop searching over all functions and instead restrict the search to the span of a finite, chosen set of basis functions $\{\chi_\mu\}_{\mu=1}^N$, usually centered on the nuclei (atomic orbital-like functions, whether Slater-type or Gaussian-type). Any molecular orbital is then written as

This is exactly the LCAO idea chemists already use qualitatively when they draw bonding and antibonding combinations of atomic orbitals — Roothaan–Hall makes it quantitative and variational. Substituting this expansion into the Fock equation and projecting onto each basis function turns a differential equation for an unknown function $\psi_i$ into a matrix equation for the unknown numbers $C_{\mu i}$. The infinite-dimensional variational search collapses to an $N$-dimensional one, and $N$-dimensional problems are exactly what matrices and computers are built to solve. Three matrices carry all the physics. The overlap matrix $S_{\mu\nu} = \langle \chi_\mu | \chi_\nu\rangle$ measures how much the basis functions overlap each other; it would be the identity matrix only for an orthonormal basis, which atom-centered functions never quite are. The density matrix $P_{\mu\nu}$ packages the occupied-orbital coefficients into a single object that directly gives the electron density, $\rho(\mathbf r) = \sum_{\mu\nu} P_{\mu\nu}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)$. The Fock matrix $F_{\mu\nu}$ is the matrix representation of $\hat f$ in the basis, built from one-electron (kinetic + nuclear attraction) integrals and two-electron Coulomb/exchange integrals weighted by $P$. Because $F$ depends on $P$, and $P$ depends on the coefficients that come from diagonalizing $F$, the whole system must be solved by iteration — the self-consistent field loop.

The starting point is the Fock eigenvalue equation for the exact (canonical) molecular orbitals,

where:

• $\hat f$ is the one-electron Fock operator, $\hat f = \hat h + \sum_j (2\hat J_j - \hat K_j)$ for closed-shell (restricted) Hartree–Fock,
• $\hat h$ is the core Hamiltonian (kinetic energy + nuclear attraction),
• $\hat J_j, \hat K_j$ are the Coulomb and exchange operators built from occupied orbital $j$,
• $\varepsilon_i$ is the orbital energy of $\psi_i$.

Reading: each electron feels the bare nuclei plus a mean-field Coulomb repulsion from every occupied orbital, corrected by an exchange term that removes the unphysical self-interaction between same-spin electrons and lowers the energy of parallel-spin pairs — this is the Pauli-driven piece. Expanding $\psi_i = \sum_\mu C_{\mu i}\chi_\mu$ and left-multiplying by $\chi_\nu^*$ before integrating gives the matrix form,

where:

• $F_{\nu\mu} = \langle \chi_\nu|\hat f|\chi_\mu\rangle$ is the Fock matrix,
• $S_{\nu\mu} = \langle \chi_\nu|\chi_\mu\rangle$ is the overlap matrix,
• $C$ is the $N\times N$ matrix of expansion coefficients (columns are orbitals),
• $\varepsilon$ is the diagonal matrix of orbital energies.

Reading: this is a generalized matrix eigenvalue problem — generalized because $S$ appears on the right, reflecting that basis functions on different centers are not orthogonal. If they were orthonormal, $S=I$ and this would be an ordinary eigenvalue equation. The Fock matrix elements are built explicitly from one- and two-electron integrals over the basis,

where:

• $H^{\text{core}}_{\mu\nu} = \langle\chi_\mu|-\tfrac12\nabla^2 - \sum_A Z_A/r_A|\chi_\nu\rangle$ is the one-electron core Hamiltonian matrix,
• $(\mu\nu|\lambda\sigma)$ is a two-electron repulsion integral in chemists' notation, $\iint \chi_\mu(1)\chi_\nu(1)\,r_{12}^{-1}\,\chi_\lambda(2)\chi_\sigma(2)\,d\mathbf r_1 d\mathbf r_2$,
• $P_{\lambda\sigma} = 2\sum_{i}^{\text{occ}} C_{\lambda i}C_{\sigma i}$ is the density matrix (factor 2 for doubly occupied closed-shell orbitals).

Reading: the Fock matrix is the core Hamiltonian dressed by a weighted sum over all basis-function pairs, with the weights supplied by the density matrix — this is precisely where the equations become nonlinear and self-consistency becomes unavoidable. Finally the total electronic energy, needed to judge convergence and to compare structures, is

where the factor $\tfrac12$ avoids double-counting the electron–electron repulsion, which the naive sum of orbital energies would count twice.

Example 1 — a minimal two-function basis, entirely by hand. Consider a toy diatomic-like system described by just two normalized real basis functions $\chi_1,\chi_2$ centered on two different atoms, with overlap $S_{12}=S_{21}=0.30$ (an off-round, physically reasonable value for moderately separated atoms) and $S_{11}=S_{22}=1$. Suppose a converged Fock matrix (in some energy unit, say hartree) has already been assembled as

We solve $FC = SC\varepsilon$ by first finding $\varepsilon$ from $\det(F-\varepsilon S)=0$:

Expanding,

The first product is $1.14 + 1.20\varepsilon + 0.95\varepsilon + \varepsilon^2 = \varepsilon^2 + 2.15\varepsilon + 1.14$. The second is $(0.45+0.30\varepsilon)^2 = 0.2025 + 0.27\varepsilon + 0.09\varepsilon^2$. Subtracting,

Applying the quadratic formula,

Since $\sqrt{0.1219}=0.3491$, the two orbital energies are

For $\varepsilon_1=-1.225$: solve $(F-\varepsilon_1 S)\mathbf c = 0$. The first row gives $(-1.20+1.225)C_1 + (-0.45+1.225(0.30))C_2 = 0$, i.e. $0.025\,C_1 + (-0.45+0.3675)C_2 = 0.025C_1 - 0.0825C_2=0$, so $C_1 = 3.30\,C_2$. Choosing $C_2=0.279$ and $C_1=0.921$ (unnormalized ratio scaled arbitrarily), we then enforce normalization $C^T S C = 1$:

With the ratio $C_1=3.30C_2$: $(3.30)^2C_2^2 + C_2^2 + 2(0.30)(3.30)C_2^2 = (10.89+1+1.98)C_2^2 = 13.87\,C_2^2=1$, so $C_2=0.2686$ and $C_1=0.886$. The lowest-energy (bonding-like) molecular orbital is

with both coefficients the same sign — a bonding combination, as expected for the lower root. Repeating the same substitution with $\varepsilon_2=-0.841$ gives an antibonding combination with $C_1$ and $C_2$ of opposite sign (the algebra is the same pattern and is left as a natural check: the second row must also be satisfied simultaneously, confirming the two roots are consistent).

Limiting case check: if $S_{12}\to 0$ (orthogonal basis), the secular equation reduces to $\det(F-\varepsilon I)=0$, an ordinary eigenvalue problem, and the mixing coefficients would simplify to the textbook two-level formula $\tan2\theta = 2F_{12}/(F_{11}-F_{22})$. This confirms the generalized problem correctly reduces to the familiar ordinary one when $S=I$.

Example 2 — energetic consequence: why the SCF loop is needed (HeH⁺-like density update). Suppose after one SCF cycle for a two-electron, two-basis-function system (one function on He, one on H, mimicking HeH⁺) the coefficients of the single doubly occupied orbital are $C_1=0.80$, $C_2=0.35$ (unnormalized for illustration, off-round on purpose). The density matrix element that feeds directly into the next Fock matrix is

If the two-electron integrals over this basis are, in some convenient units, $(11|11)=0.90$, $(22|22)=0.55$, $(11|22)=0.40$, $(12|12)=0.15$ (all asymmetric, chosen to avoid accidental cancellation), the Coulomb contribution to $F_{11}$ from this density is

Taking $(11|12)=(11|21)=0.20$ for concreteness, this evaluates to

This number changes $F_{11}$, which changes $\varepsilon$ and $C$ on the next diagonalization, which changes $P_{11}$ again — precisely the nonlinear feedback loop that the derivation above predicted. Only when $P_{11}$ (and every other density matrix element) stabilizes between cycles has self-consistency been reached, and only then is the energy formula $E_{\text{elec}}=\tfrac12\sum P_{\mu\nu}(H^{\text{core}}_{\mu\nu}+F_{\mu\nu})$ meaningful as the Hartree–Fock energy.

Example 3 — orthogonalizing a non-orthogonal basis (Löwdin transformation), exam-style. A common exam question asks students to convert $FC=SC\varepsilon$ into an ordinary eigenvalue problem. Given the same overlap matrix as Example 1, $S=\begin{pmatrix}1 & 0.30\\0.30 & 1\end{pmatrix}$, find $S^{-1/2}$ and the transformed Fock matrix $F'$. First diagonalize $S$: its eigenvalues solve $\det(S-\lambda I)=0 \Rightarrow (1-\lambda)^2 - 0.09=0 \Rightarrow \lambda = 1\pm0.30$, giving $\lambda_1=1.30$ and $\lambda_2=0.70$, with (by symmetry of this particular symmetric $2\times2$ matrix) eigenvectors $\tfrac{1}{\sqrt2}(1,1)$ and $\tfrac1{\sqrt2}(1,-1)$. Then

with $\lambda_1^{-1/2}=1/\sqrt{1.30}=0.8771$ and $\lambda_2^{-1/2}=1/\sqrt{0.70}=1.1952$. Carrying out the matrix product,

So $S^{-1/2} = \begin{pmatrix}1.0362 & -0.1591\\ -0.1591 & 1.0362\end{pmatrix}$. Using the Fock matrix from Example 1, the transformed (orthogonalized) Fock matrix is $F' = S^{-1/2}FS^{-1/2}$. Computing the first product, $FS^{-1/2}$:

By the analogous computation for the second row ($F_{21}=F_{12}=-0.45$, $F_{22}=-0.95$):

Then $F' = S^{-1/2}(FS^{-1/2})$, and its diagonal entries are

Diagonalizing this ordinary (symmetric, orthogonal-basis) $F'$ by the standard rotation formula reproduces the same physical orbital energies found in Example 1 ($\varepsilon_1=-1.225$, $\varepsilon_2=-0.841$ hartree, matching up to the rounding carried through this hand calculation), because Löwdin orthogonalization is an exact similarity-type transformation, not an approximation — it only changes the basis representation, never the physics.

George G. Hall, working independently in Cambridge and publishing essentially the same matrix equations in the same year, is honored by attaching his name alongside Roothaan's; the equations are properly a joint discovery even though they were derived without direct collaboration. The near-simultaneous, independent arrival at identical matrix equations is itself informative: once the LCAO ansatz for The Self-Consistent-Field Procedure was combined with standard linear algebra, the matrix formulation was, in a real sense, inevitable — the physics and the available mathematics of the era pointed to exactly one destination.