SCF convergence theory: DIIS, SOSCF, and augmented-Hessian methods

The SCF chapter showed the convergence controls as a practical menu. This chapter derives them: the unified picture that makes every strategy a special case of “drive the orbital gradient to zero,” the full CDIIS/EDIIS/ADIIS construction, the stabilizers (damping, level shift, Fermi smearing), and the second-order family (SOSCF, the augmented Hessian, TRAH) that finishes what first-order methods cannot.

The unifying picture: SCF as unconstrained optimization

Parametrize any new set of orbitals reachable from the current \(\mathbf C\) by a unitary rotation,

\[ \mathbf C(\kappa) = \mathbf C\,\exp(\kappa), \qquad \kappa_{ai}=-\kappa_{ia}^{*},\quad \kappa_{ij}=\kappa_{ab}=0, \]

restricted to the occupied–virtual block (\(i,j\) occupied, \(a,b\) virtual) because rotations within the occupied or virtual subspace leave the Slater determinant, and hence the energy, invariant — only the occ–virt block is non-redundant. Expand the energy in \(\kappa\) to second order:

\[ E(\kappa) = E_0 + \mathbf g^{\mathsf T}\kappa + \tfrac12\,\kappa^{\mathsf T}\mathbf H\,\kappa + \mathcal O(\kappa^3). \]

For RHF the orbital gradient and orbital Hessian are

\[ g_{ai} = 4F_{ai}, \qquad H_{ai,bj} = \delta_{ij}F_{ab} - \delta_{ab}F_{ij} + 4(ai|bj) - (ab|ij) - (aj|bi), \]

with \(F_{ai}\) the occ–virt block of the MO Fock matrix and \((pq|rs)\) the MO two-electron integrals (UHF sums per spin channel; ROHF needs the three-block closed/open/virtual generalization, with the full cross-block Fock commutators worked out so the Hessian remains exact). The stationarity condition \(\mathbf g=0\) is the SCF convergence condition, and it has an exactly equivalent AO-basis form: the commutator \(\mathbf{FDS}-\mathbf{SDF}\) vanishes at self-consistency, because that commutator is precisely the occ–virt gradient block written in the AO representation rather than the MO one. This single equivalence is why every method below — however different its bookkeeping — is really “a different approximation to a step that drives \(\mathbf g\to0\)”:

Method family	Uses	Character
CDIIS / EDIIS / ADIIS	\(\mathbf g\) only (extrapolated across cycles)	first-order
Damping / level shift / smearing	\(\mathbf g\), stabilized	first-order + regularization
SOSCF	\(\mathbf g\) + an approximate \(\mathbf H^{-1}\) (BFGS)	quasi-second-order
Augmented Hessian (QC-SCF) / TRAH	\(\mathbf g\) + the exact \(\mathbf H\), via Hessian–vector products	second-order

No method ever builds \(\mathbf H\) explicitly (\(\mathcal O(n_{\text{ao}}^4)\) storage) — every second-order strategy evaluates only the Hessian–vector product \(\sigma=\mathbf H\kappa\), at the cost of one Fock-like contraction against a transition density rather than the true one. This response engine — exact for RHF/UHF/ROHF, and including the exact KS exchange–correlation kernel \(f_{xc}\) and the PCM reaction-field response when active — is the single piece of machinery every second-order method and the analytic Hessian (next chapters) build on.

First-order methods: extrapolating the gradient

CDIIS (Pulay)

Direct inversion of the iterative subspace (Pulay, 1980/1982) extrapolates a better Fock matrix from the last \(n\) cycles, using each cycle’s commutator as an error signal.

Definition 3 (The DIIS error vector)

In an orthonormal basis (\(\mathbf X=\mathbf S^{-1/2}\)), the error at cycle \(i\) is $\( \mathbf e_i = \mathbf X^{\mathsf T}\big(\mathbf F_i\mathbf D_i\mathbf S - \mathbf S\mathbf D_i\mathbf F_i\big)\mathbf X, \)$ which vanishes exactly at self-consistency (it is the AO-basis orbital gradient).

Algorithm 1 (CDIIS extrapolation)

Input: the last \(n\) Fock/density pairs \(\{\mathbf F_i,\mathbf D_i\}\). Output: an extrapolated \(\bar{\mathbf F}\) to diagonalize next.

1. Build the error \(\mathbf e_i\) for each history slot (Definition above).

2. Form the Gram matrix \(B_{ij}=\operatorname{Tr}(\mathbf e_i^{\mathsf T}\mathbf e_j)\).

3. Minimize \(\|\sum_i c_i\mathbf e_i\|^2\) subject to \(\sum_i c_i=1\). Introducing a Lagrange multiplier \(\lambda\) for the constraint and setting the gradient of the Lagrangian to zero gives the bordered linear system $$

   (1)\[\begin{pmatrix}\mathbf B & -\mathbf 1\\ -\mathbf 1^{\mathsf T} & 0\end{pmatrix}\]

   \begin{pmatrix}\mathbf c\ \lambda\end{pmatrix} = \begin{pmatrix}\mathbf 0\ -1\end{pmatrix}. $$

4. Extrapolate \(\bar{\mathbf F}=\sum_i c_i\mathbf F_i\) and diagonalize it in place of the latest \(\mathbf F_n\).

CDIIS converges superlinearly near the solution — its whole strength — but minimizing an error norm is not the same as lowering the energy: nothing in the construction prevents an extrapolated \(\bar{\mathbf F}\) from a step that raises the energy when the history is still far from self-consistent. In practice this shows up as instability or outright divergence starting from a poor guess, which motivates the energy-based alternatives below. Two further practical safeguards matter: an ill-conditioned \(\mathbf B\) (near-degenerate error vectors) needs the oldest history vector dropped or the diagonal regularized, and the coefficients \(|c_i|\) need an upper bound — an unconstrained least-squares solution can occasionally produce wild extrapolation coefficients that overshoot badly.

EDIIS and ADIIS: constrained energy minimization

EDIIS (Kudin, Scuseria, Cancès, 2002) instead minimizes an energy model directly. For a convex combination of densities \(\mathbf D(\mathbf c)=\sum_i c_i\mathbf D_i\) with \(c_i\ge0,\sum_i c_i=1\), and because the HF energy is a quadratic functional of the density, the trial energy has the closed form

\[ E^{\text{EDIIS}}(\mathbf c) = \sum_i c_i E_i - \tfrac12\sum_{i,j}c_ic_j\langle\mathbf F_i-\mathbf F_j,\, \mathbf D_i-\mathbf D_j\rangle, \]

with \(E_i\) the total energy at cycle \(i\) and \(\langle\cdot,\cdot\rangle\) the Frobenius inner product. This is a small (history-size), box-constrained quadratic program — cheap to solve exactly with an active-set or projected-gradient method — and the constraint \(c_i\ge0\) is exactly what makes EDIIS robust far from convergence: a negative coefficient is precisely what would let CDIIS overshoot, and EDIIS structurally forbids it.

ADIIS (Hu, Yang, 2010) is EDIIS’s DFT-safe refinement: rather than assume a quadratic energy (which KS-DFT is not, because \(E_{xc}\) is not quadratic in the density), it builds an augmented Roothaan–Hall model referenced to the most recent iterate \(n\),

\[ E^{\text{ADIIS}}(\mathbf c) = 2\sum_i c_i\langle\mathbf D_i-\mathbf D_n,\,\mathbf F_n\rangle + \sum_{i,j}c_ic_j\langle\mathbf D_i-\mathbf D_n,\,\mathbf F_j-\mathbf F_n\rangle, \]

again minimized under \(c_i\ge0,\sum c_i=1\). ADIIS and EDIIS coincide exactly for Hartree–Fock (where the quadratic assumption is exact) and ADIIS is the more robust choice for DFT. qc-rs’s algorithm="auto" blends these energy-based extrapolations with CDIIS by monitoring the error norm \(\|\mathbf e_n\|\): energy-based extrapolation while \(\|\mathbf e_n\|\) is large, a linear interpolation through the mid-range, and pure CDIIS once \(\|\mathbf e_n\|\) drops below a tight threshold, where CDIIS’s superlinear convergence is unmatched.

Stabilizers

Three further tools steady any of the strategies above without replacing them.

Damping linearly mixes the new and previous density, $\( \tilde{\mathbf D}^{(n)} = (1-\alpha)\,\mathbf D^{(n)}_{\text{new}} + \alpha\,\mathbf D^{(n-1)},\qquad 0\le\alpha<1, \)\( a blunt but effective brake on early oscillation (\)\alpha\approx0.7\( is a typical static value, released once the error drops below a threshold; a dynamic Hehenberger–Zerner scheme adjusts \)\alpha$ from the energy/gradient ratio each cycle).

Level shifting (Saunders–Hillier, 1973) raises the virtual-orbital energies by a constant \(b\) to suppress spurious occupied–virtual mixing near a small HOMO–LUMO gap. In the MO basis this is simply \(\tilde F_{pq}=F_{pq}+b\,\delta_{pq}\) for \(p,q\in\text{virtual}\); the equivalent AO-basis form, using the virtual-space projector \(\mathbf Q=\mathbf S^{-1}-\mathbf D\), $\( \tilde{\mathbf F} = \mathbf F + b\,\mathbf S\mathbf Q\mathbf S = \mathbf F + b\,(\mathbf S - \mathbf S\mathbf D\mathbf S), \)$ lets it be added directly to the AO Fock matrix. A shift too large slows the late stage of convergence (it artificially inflates the gap it is meant to only stabilize), so it is released as the error shrinks.

Fermi smearing replaces the integer Aufbau occupation with a finite-temperature Fermi–Dirac distribution, $\( f_i = \Big[1+\exp\big((\varepsilon_i-\mu)/k_BT_{\text{el}}\big)\Big]^{-1}, \)\( solving for the chemical potential \)\mu\( by bisection/Newton so that \)\sum_i w_if_i=N_{\text{el}}\( (\)f_i\( is monotonic in \)\mu\(, so the root is unique). The density becomes \)\mathbf D=\sum_if_i|\mathbf c_i\rangle\langle \mathbf c_i|\(, and because fractional occupation is a genuine free-energy minimization, convergence should track the electronic free energy \)\( A = E - T_{\text{el}}S_{\text{el}}, \qquad S_{\text{el}} = -k_B\sum_i\big[f_i\ln f_i+(1-f_i)\ln(1-f_i)\big], \)\( not the bare energy. Smearing is exactly the fix for **near-degenerate or metallic** systems where integer occupation flips orbitals in and out every cycle, causing persistent oscillation; a common annealing schedule starts hot and cools \)T_{\text{el}}\to0$ as the SCF converges, recovering an ordinary integer-occupation solution at the end.

Second-order methods: using the exact Hessian

When first-order methods trail (typically once \(\|\mathbf e\|\sim10^{-3}\)–\(10^{-4}\)) or fail outright, the second-order family takes an explicit Newton-flavoured step using \(\mathbf H\).

SOSCF (Chaban–Schmidt–Gordon, 1997) never forms \(\mathbf H\) at all — it approximates its inverse with BFGS quasi-Newton updates from an (s,y)-pair history, seeded by the cheap diagonal guess \(H^{(0)}_{ai,ai}\approx4(\varepsilon_a-\varepsilon_i)\), and steps \(\kappa=-\mathbf H^{-1}\mathbf g\). With no Hessian–vector product at all it is the cheapest second-order option, converging superlinearly near the solution but unstable started far away — the reason it is a DIIS finisher, engaged once the gradient is already small, rather than a starting strategy.

The augmented Hessian / QC-SCF (Bacskay, 1981) confronts Newton’s real failure mode directly: a plain step \(\mathbf H\kappa=-\mathbf g\) does not even descend if \(\mathbf H\) has a negative eigenvalue, which is common far from convergence (it is exactly the signature tested for in stability analysis). The fix is to solve the augmented eigenvalue problem for its lowest eigenpair, $$

(2)\[\begin{pmatrix}0 & \alpha\mathbf g^{\mathsf T}\\ \alpha\mathbf g & \mathbf H\end{pmatrix}\]

(3)\[\begin{pmatrix}1\\ \tilde\kappa\end{pmatrix} = \mu\begin{pmatrix}1\\ \tilde\kappa\end{pmatrix}\]

;\Longrightarrow; (\mathbf H-\mu\mathbf I)\kappa=-\mathbf g, $\( whose eigenvalue \)\mu<0\( is an **automatic level shift**: \)\mathbf H-\mu\mathbf I\( is positive-definite by construction, so the resulting step descends *even when \)\mathbf H\( itself is not*. Because only the lowest eigenpair is needed, a **Davidson** iteration finds it by repeatedly evaluating the augmented-matrix product on the fly — never assembling \)\mathbf H$.

TRAH (trust-region augmented Hessian; Helmich-Paris, 2021) adds an explicit trust region to the same augmented problem, rescaling \(\mathbf H\) by a trial step length \(\lambda\ge1\), $$

(4)\[\begin{pmatrix}0 & \mathbf g^{\mathsf T}\\ \mathbf g & \mathbf H/\lambda\end{pmatrix}\]

(5)\[\begin{pmatrix}1\\ x\end{pmatrix} = \varepsilon\begin{pmatrix}1\\ x\end{pmatrix}\]

;\Longrightarrow; \kappa=\frac{x}{\lambda} = -(\mathbf H-\lambda\varepsilon\mathbf I)^{-1}\mathbf g, $\( choosing \)\lambda\( so the resulting \)|\kappa|\( matches a target trust radius, then growing or shrinking that radius by a \)\rho\(-test comparing the *predicted* energy drop (from the model) against the *actual* one. A key implementation identity — descent is guaranteed by \)\mathbf g^{\mathsf T}\kappa=\varepsilon/\lambda\le0\( without needing \)\mathbf H\( to be positive-definite — lets the trust radius adapt without recomputing the expensive Hessian–vector products at each trial \)\lambda$. TRAH’s adaptive radius is what makes it reliable on the hardest cases (an orbitally near-degenerate ROHF radical, say) where a fixed step can overshoot into the wrong electronic basin.

Convergence criteria

Production SCF codes require several criteria simultaneously, not just the energy:

Criterion	Symbol	Typical tight threshold
Energy change	$	\Delta E
Density RMS change	\(\operatorname{RMS}(\Delta\mathbf D)\)	\(5\times10^{-9}\)
Density max change	$\max	\Delta\mathbf D
DIIS error (commutator)	\(|\mathbf e|\)	\(5\times10^{-7}\)
Orbital gradient	\(|\mathbf g|\)	\(10^{-5}\)

A second-order method is naturally judged on \(\|\mathbf g\|\) alone (it is the quantity the method drives to zero directly); first-order methods benefit from the fuller multi-criterion AND, since an energy plateau alone can mask an SCF that has not actually reached self-consistency in the density.

The recommended strategy ladder

        flowchart TD
    G["Guess (SAD / SAP)"] --> F["Far field: ||e|| > 0.1"]
    F -->|"EDIIS (HF) / ADIIS (DFT)<br/>+ damping/level-shift if needed"| M["Mid field: 1e-3 < ||e|| <= 0.1"]
    M -->|"blend of energy-model + CDIIS"| N["Near field: ||e|| <= 1e-3"]
    N -->|"pure CDIIS (fastest locally)"| T{"trailing?"}
    T -->|yes| S["SOSCF (superlinear finisher)"]
    T -->|no, converged| D["Done"]
    N -->|"stall / oscillation"| SM["Fermi smearing + anneal T -> 0"]
    SM -->|"still hard"| AH["Augmented Hessian / TRAH<br/>(globally convergent)"]
    S --> D
    AH --> D
    

Every rung is “a different approximation to \(\mathbf g\to0\),” which is why the ladder always converges to the same energy regardless of which rungs a given molecule needs — exactly the empirical demonstration in the SCF chapter. algorithm="xqc" automates the common case (plain DIIS for a bounded number of cycles, falling back to the augmented-Hessian ladder only if unconverged); algorithm="yqc" is the same idea tuned for orbitally-degenerate open shells (ROHF/ROKS), which stall a pure first-order DIIS long before reaching the near-field regime and so default to this escalation rather than plain diis.

Exercise 4

1. Why must EDIIS’s coefficients satisfy \(c_i\ge0\) while CDIIS’s do not, and how does that constraint alone explain EDIIS’s robustness far from convergence?

2. The augmented-Hessian eigenvalue equation guarantees descent even when \(\mathbf H\) has a negative eigenvalue. Identify the one quantity in the construction that makes this true, and explain in one sentence why.

3. TRAH rescales \(\mathbf H\) by \(1/\lambda\) rather than recomputing the Hessian–vector product at each trial trust radius. Why does this matter for TRAH’s cost per macro-iteration?

Solution to Exercise 4

1. \(c_i\ge0\) restricts \(\mathbf D(\mathbf c)=\sum_ic_i\mathbf D_i\) to a convex combination of the history densities — a point inside the simplex spanned by them. CDIIS’s unconstrained coefficients can be negative, which is algebraically an extrapolation beyond the historical densities — exactly the kind of overshoot that destabilizes CDIIS far from convergence, and exactly what \(c_i\ge0\) forbids.

2. The automatic level shift \(\mu\) (the augmented problem’s lowest eigenvalue). Because \((\mathbf H-\mu\mathbf I)\kappa=-\mathbf g\) and \(\mu\) is by construction at or below \(\mathbf H\)’s lowest eigenvalue, \(\mathbf H-\mu\mathbf I\) is positive-semidefinite regardless of \(\mathbf H\)’s own sign structure, guaranteeing a descent direction.

3. Hessian–vector products (via the response engine) are the expensive part of a macro-iteration. By keeping the Krylov/Davidson subspace matrices split as \(\tilde M = \tilde M_1 + \tilde M_2/\lambda\), TRAH can retry several trial trust radii \(\lambda\) using only cheap linear algebra on the already-computed subspace matrices — no new \(\sigma=\mathbf H\kappa\) evaluation per trial.

The response engine behind \(\sigma=\mathbf H\kappa\) reappears twice more in this manual: as the coupled-perturbed equations behind the analytic Hessian, and as the operator whose eigenvalues stability analysis tests directly.