# SCF convergence theory: DIIS, SOSCF, and augmented-Hessian methods
The [SCF chapter](../20-guide/scf.md) 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.
:::{prf:definition} The DIIS error vector
:label: def-diis-error
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).
:::
:::{prf:algorithm} CDIIS extrapolation
:label: alg-cdiis-full
**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
$$
\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](../20-guide/scf.md)). The fix is to solve the **augmented eigenvalue problem** for its
lowest eigenpair,
$$
\begin{pmatrix}0 & \alpha\mathbf g^{\mathsf T}\\ \alpha\mathbf g & \mathbf H\end{pmatrix}
\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$,
$$
\begin{pmatrix}0 & \mathbf g^{\mathsf T}\\ \mathbf g & \mathbf H/\lambda\end{pmatrix}
\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|$ | $10^{-8}$ |
| Density RMS change | $\operatorname{RMS}(\Delta\mathbf D)$ | $5\times10^{-9}$ |
| Density max change | $\max|\Delta\mathbf D|$ | $10^{-7}$ |
| 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
```{mermaid}
flowchart TD
G["Guess (SAD / SAP)"] --> F["Far field: ||e|| > 0.1"]
F -->|"EDIIS (HF) / ADIIS (DFT)
+ 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
(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](../20-guide/scf.md)
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}
:label: ex-scf-theory
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} ex-scf-theory
:class: dropdown
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](analytic-hessian-thermo.md), and as the operator whose eigenvalues
[stability analysis](../20-guide/scf.md) tests directly.