# SCF: Hartree–Fock & KS-DFT
You have a molecule, a basis, and a starting guess. This chapter runs the calculation that turns them into
orbitals and an energy: the **self-consistent field (SCF)**, which solves both **Hartree–Fock** and
**Kohn–Sham DFT**. It is the workhorse every later method builds on, so this is the longest chapter in the
guide — but the controls are few, and most runs need only two of them (`ref` and `xc`).
## Theory recap: what the SCF solves
From [Part II](../10-foundations/hartree-fock.md): both HF and KS-DFT reduce to a one-electron eigenvalue
problem in the finite basis, the **Roothaan equations**
$$
\mathbf{F}\,\mathbf{C} = \mathbf{S}\,\mathbf{C}\,\boldsymbol\varepsilon ,
$$
where the Fock (or Kohn–Sham) matrix $\mathbf F$ depends on the density $\mathbf D$ built from the very
orbitals $\mathbf C$ it produces. Because of that circularity the equations are solved **iteratively** —
build $\mathbf F$, diagonalize, rebuild the density, repeat until it stops changing ([the initial-guess
chapter](initial-guess.md) started that loop). qc-rs handles the whole cycle; your job is to specify **which
electronic structure** to solve for (the reference and, for DFT, the functional) and, occasionally, **how to
steer a hard convergence**.
## The basic run
```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
done.scf.energy # -76.026794 total energy (hartree)
done.scf.converged # True
done.scf.ncycle # 9 SCF iterations used
done.scf.energy_elec # -85.220707 electronic energy (= total − nuclear repulsion)
done.scf.energy_components # {'core': -123.1493, 'coulomb': 46.9052, 'exchange': -8.9766}
```
`.scf(...)` adds a pending SCF step; `.run()` executes it (the functional `qc.scf(mychk, ...)` is
equivalent). The **`energy_components`** break the electronic energy into its physical pieces — for HF the
one-electron `core` (kinetic + nuclear attraction), the `coulomb` repulsion $J$, and the `exchange` $K$; for
DFT the last becomes an `xc` term (below).
## Choosing the reference: RHF, UHF, ROHF
The reference is set by **`ref`** together with the spin multiplicity **`spin` = 2S+1** (from
[`qc.chk.new`](molecular-input.md)). The three references differ in how they constrain the spin-orbitals
([Part II](../10-foundations/hartree-fock.md)):
| `ref` | closed shell (`spin=1`) | open shell (`spin>1`) |
|---|---|---|
| `"auto"` *(default)* | **RHF** / RKS | **UHF** / UKS |
| `"r"` | **RHF** / RKS | **ROHF** / ROKS |
| `"u"` | UHF / UKS | UHF / UKS |
| `"ro"` | **error** | **ROHF** / ROKS |
**`ref="auto"` (the default)** picks *restricted* for a closed shell and *unrestricted* for an open shell —
so a radical or triplet runs UHF without your asking (UHF is the conventional open-shell default). `ref="ro"`
is the explicit request for restricted open-shell; it **errors on a closed shell** rather than silently
running RHF, as a guard against a wrong `charge`/`spin`.
### UHF vs ROHF: a trade-off
For an open shell you choose between two references, and they give *different* energies:
```python
ch3 = "C 0 0 0; H 0 1.079 0; H 0.934 -0.539 0; H -0.934 -0.539 0" # methyl radical, a doublet
u = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u").run().scf
ro = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="ro").run().scf
print(u.energy, ro.energy) # -39.563802 -39.559636
```
- **UHF** lets the α and β electrons occupy *different* spatial orbitals. That extra freedom gives a **lower
energy** (−39.563802 vs −39.559636), but the wavefunction is no longer a pure spin state — it suffers
**spin contamination**, measurable as $\langle S^2\rangle$ drifting above the exact value ($S(S+1)=0.75$
for a doublet):
```python
qc.prop.spin.s_squared(qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u").run())
# 0.7612 (exact doublet = 0.75; the small excess is the contamination)
```
- **ROHF** forces paired electrons to share a spatial orbital, so it stays **spin-pure** ($\langle
S^2\rangle$ exact) at the cost of a slightly higher energy.
Use **UHF/UKS** as the default for radicals, and **ROHF/ROKS** when a spin-pure reference matters (e.g. as a
basis for later correlation, or to avoid contamination artefacts).
:::{important} Near-degenerate open shells → UHF, not ROHF
A near-degenerate or delocalized-hole state — e.g. a symmetric cluster cation where the hole spreads over
equivalent sites — is effectively multireference and **cannot** be represented by a single-determinant ROHF;
it simply will not converge. UHF reaches it by symmetry-breaking (localizing the hole) and is the right
choice. ROHF/ROKS is for well-defined high-spin or doublet radicals.
:::
## Hartree–Fock or KS-DFT: `xc=`
The *same* SCF machinery runs both methods; **`xc=`** is the only switch:
- **`xc=None`** (the default) or **`xc="hf"`** → **Hartree–Fock**.
- **`xc="pbe"`**, **`"b3lyp"`**, … → **Kohn–Sham DFT** with that libxc functional (LDA / GGA / meta-GGA /
hybrids; the theory and Jacob's ladder are in [Part II](../10-foundations/dft-kohn-sham.md)).
```python
pbe = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe").run().scf
print(pbe.energy) # -76.333409
print(dict(pbe.energy_components)) # {'core': -123.1775, 'coulomb': 46.9303, 'xc': -9.2801}
```
Notice the components now report an **`xc`** term (exchange–correlation) instead of HF's `exchange` — the one
structural difference between an HF and a KS run.
### The DFT grid
KS-DFT evaluates the exchange–correlation energy by **numerical quadrature** on a molecular grid, chosen with
**`grid=`** (ORCA-DefGrid-like pruned tiers):
| `grid=` | radial × angular | use |
|---|---|---|
| `coarse` | 50 × 194 | draft / quick |
| `medium` *(default)* | 60 × 302 | production (≈ ORCA DefGrid2) |
| `fine` | 75 × 434 | tight |
| `ultrafine` | 200 × 5810 | reference / max accuracy |
The default `medium` is calibrated within ≲0.4 µEh/atom of the dense `ultrafine` reference, so you rarely
need to change it; step up to `fine` for tight energy differences or when a meta-GGA looks grid-sensitive.
(HF has no grid — it is analytic.)
## Theory: the SCF as an optimization problem
Before the convergence *strategies* make sense, it helps to see what they are all approximating. Rewrite the
Roothaan equations as an **unconstrained optimization** over orbital rotations. Any new set of orbitals
$\mathbf C(\kappa)$ reachable from the current $\mathbf C$ by a unitary rotation can be written
$$
\mathbf C(\kappa) = \mathbf C\,\exp(\kappa), \qquad \kappa_{ai} = -\kappa_{ia}^{*},\quad \kappa_{ii}=\kappa_{aa}=0,
$$
where $\kappa$ is an anti-Hermitian matrix restricted to the **occupied–virtual** block (indices $i,j$
occupied; $a,b$ virtual) — rotations *within* the occupied or virtual space leave the Slater determinant, and
hence the energy, unchanged, so only the occ–virt block is non-redundant. Expanding 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),
$$
gives an **orbital gradient** $\mathbf g$ and an **orbital Hessian** $\mathbf H$. For RHF,
$$
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. The
gradient vanishing, $\mathbf g=0$, is exactly the SCF convergence condition — and it is equivalent to the
familiar AO commutator $\mathbf{FDS}-\mathbf{SDF}=0$, since that commutator *is* the occ–virt block of the
gradient written in the AO basis. This single fact is the thread tying every strategy below together:
:::{important} One picture, many algorithms
**Every convergence strategy is a different approximation to "take a step that drives $\mathbf g\to 0$."**
First-order methods (DIIS) use only $\mathbf g$ (extrapolated across cycles); second-order methods (SOSCF, the
augmented-Hessian family) also use $\mathbf H$, exactly or approximately, to take a genuine Newton-like step.
Building $\mathbf H$ explicitly costs $\mathcal O(n_{\text{ao}}^4)$ and is never done — instead every
second-order method evaluates only the **Hessian–vector product** $\sigma=\mathbf H\kappa$, which costs the
same as one extra Fock build (a contraction with a "transition density" rather than the true one).
:::
## Theory: DIIS — direct inversion of the iterative subspace
**DIIS** (Pulay, 1980/1982) is the default first-order accelerator. Instead of taking the current Fock
matrix at face value, it **extrapolates** a better one from the last several cycles' worth of Fock matrices,
using their commutator errors as a guide.
:::{prf:algorithm} CDIIS (Pulay) step
:label: alg-cdiis
**Input:** the last $n$ Fock/density pairs $\{F_i, D_i\}$; the orthogonalizer $X=\mathbf S^{-1/2}$.
**Output:** an extrapolated Fock matrix $\bar F$ to diagonalize.
1. For each history slot $i$, form the **error matrix** in an orthonormal basis,
$$
e_i = X^{\mathsf T}\big(F_i D_i S - S D_i F_i\big) X .
$$
$e_i \to 0$ at self-consistency (it is the AO form of the orbital gradient above).
2. Build the $n\times n$ Gram matrix $B_{ij} = \operatorname{Tr}\!\big(e_i^{\mathsf T} e_j\big)$.
3. Solve the Lagrange-bordered linear system for the extrapolation coefficients $\mathbf c$ (constrained to
$\sum_i c_i = 1$, enforced by the multiplier $\lambda$):
$$
\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 F = \sum_i c_i F_i$ and diagonalize *that* instead of the latest $F_n$.
:::
CDIIS minimizes the error norm $\|\sum_i c_i e_i\|$, but that norm is not the energy — nothing stops the
extrapolated $\bar F$ from a step that *raises* the energy far from convergence, which is why CDIIS alone can
be unstable early on. The **`auto`** strategy therefore blends in an **energy-based** extrapolation
(EDIIS/ADIIS — minimizing a quadratic model of the energy itself, constrained to $c_i\ge0$, so it cannot
overshoot) while the error is large, and hands off to pure CDIIS once $\|e\|$ drops below a threshold, where
CDIIS's superlinear local convergence is fastest. This is the "far-field energy-robust, near-field
error-robust" handoff described in the [convergence-theory chapter](../10-foundations/scf-convergence-theory.md)
of Part II.
## Theory: second-order convergence — SOSCF, the augmented Hessian, and TRAH
When DIIS trails (typically once $\|e\|\sim10^{-3}$–$10^{-4}$) or fails to converge at all, the second-order
family takes an explicit step in $\kappa$ using the Hessian:
- **SOSCF** (Chaban–Schmidt–Gordon) never builds $\mathbf H$ at all — it approximates its *inverse* with
**BFGS** quasi-Newton updates, starting from the cheap diagonal guess
$H^{(0)}_{ai,ai}\approx 4(\varepsilon_a-\varepsilon_i)$, and steps $\kappa=-\mathbf H^{-1}\mathbf g$. It is
the cheapest second-order option (no Hessian–vector products at all) and converges superlinearly once
close to the solution, but can be unstable far away — hence `soscf` is best as a DIIS *finisher*, not a
starting strategy.
- **QC-SCF / the augmented Hessian** (Bacskay) fixes Newton's failure mode directly: a plain Newton step
$\mathbf H\kappa=-\mathbf g$ does not even descend if $\mathbf H$ has a negative eigenvalue (common far from
convergence). Instead, solve the **augmented eigenvalue problem**
$$
\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,
$$
for its **lowest** eigenpair. The eigenvalue $\mu<0$ acts as an automatic level shift, so
$\mathbf H-\mu\mathbf I$ is always positive-definite and the step always descends — even when $\mathbf H$
itself is not. Only the lowest eigenpair is needed, found by a **Davidson** iteration that evaluates the
augmented-matrix product on the fly (never building $\mathbf H$).
- **TRAH** (trust-region augmented Hessian) adds an explicit **trust radius**: the same augmented problem is
solved with $\mathbf H$ rescaled by a trial step length $\lambda\ge1$, chosen so the resulting $\|\kappa\|$
matches a target radius — enlarging the radius when a step is accepted (`ρ`-test, comparing predicted vs.
actual energy drop) and shrinking it otherwise. This is what makes `trah` reliable on the hardest
open-shell cases (an orbitally near-degenerate ROHF radical, say) where a fixed step can overshoot into the
wrong electronic state.
```{mermaid}
flowchart TD
G["Guess density D₀ (sad, gwh, ...)"] --> F["Build Fock F"]
F --> E["Error e = FDS − SDF"]
E --> Q{"‖e‖ vs threshold"}
Q -->|"large (far)"| ED["EDIIS / ADIIS
energy-model step"]
Q -->|"small (near)"| CD["CDIIS
Pulay extrapolation"]
Q -->|"stalling"| SO["2nd-order finisher
SOSCF / QC-SCF / TRAH"]
ED --> F
CD --> F
SO --> F
Q -->|"converged"| D["Done: E, C, D"]
```
This is exactly the ladder `algorithm="auto"` walks for RHF/UHF/RKS/UKS (plain DIIS is usually enough near
equilibrium, so the energy-model branch and the 2nd-order branch rarely engage); ROHF/ROKS default to `yqc`
because an orbitally-degenerate open shell stalls a pure first-order DIIS long before reaching the near-field
regime. The full derivation — including why the Hessian–vector product never materializes $\mathbf H$, the
EDIIS/ADIIS energy functionals, and the ROHF three-block Hessian — is in
[SCF convergence theory](../10-foundations/scf-convergence-theory.md).
## Convergence: steering the SCF
With the theory above in hand, the practical controls are two dials: a **strategy** selector `algorithm=`
(*which* of the pictures above to use) and a few **stabilizer** kwargs that steady any of them. Like the
[initial guess](initial-guess.md), **the convergence controls change only the *path*, never the converged
answer** — every setting below reaches the same energy.
The default, `algorithm="auto"`, walks the ladder in the diagram above from the SAD guess, which converges
most molecules near equilibrium in well under 20 cycles. You usually change nothing. When a run *is* hard,
pick by **symptom**:
| Symptom | Reach for | How |
|---|---|---|
| Hard / oscillating — want robust 2nd order | QC-SCF (augmented Hessian) | `algorithm="qc"` |
| Same, with an adaptive trust region | TRAH | `algorithm="trah"` |
| Robust DIIS without thinking | XQC (DIIS + a 2nd-order safety net) | `algorithm="xqc"` |
| Only the tail is slow | SOSCF | `algorithm="soscf"` |
| Oscillating far from convergence | damping | `damping=0.6` |
| Near-degenerate / small-gap / metallic | Fermi smearing | `smearing=0.01` |
| HOMO/LUMO swapping each cycle | level shift | `level_shift=0.3` |
Seeing the theory in action: the augmented-Hessian family reaches the same energy in far fewer cycles on a
hard open-shell case, at the cost of a Hessian–vector product per cycle instead of a plain Fock build:
```python
for algo in ("diis", "qc", "trah"):
c = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u", algorithm=algo).run().scf
print(algo, c.ncycle, round(c.energy, 6))
# diis 11 -39.563802
# qc 6 -39.563802
# trah 6 -39.563802
```
Same energy, but the augmented-Hessian methods reach it in 6 cycles instead of 11.
### Thresholds
The SCF stops when **both** the energy change and the orbital-gradient RMS fall below tolerance:
| Threshold | default | set via |
|---|---|---|
| energy `ΔE` | `1e-9` Ha | `scf(conv_tol=...)` or `conv_preset=` |
| gradient RMS of the `[F, DS]` commutator | `1e-6` | `iop={"scf.conv_tol_grad": ...}` |
Use **`conv_preset="tight"`** / `"loose"` for a quick tighten/loosen, or `conv_tol=` for the energy directly.
The IOP keys behind the fine-grained knobs are in the [reference](../40-reference/iop.md).
## Reading and diagnosing a run
After `.run()`, the `scf` accessor exposes the result; **`run(log=...)`** (from the
[quickstart](../00-intro/quickstart.md)) streams the live cycle-by-cycle transcript that lets you *watch*
convergence:
```python
done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(log="stdout")
done.scf.energy, done.scf.converged, done.scf.ncycle # the headline numbers
dict(done.scf.energy_components) # where the energy comes from
done.log() # replay the stored transcript without recomputing
done.show("result") # a rendered snapshot of state + results
```
If `converged` is `False`, look at the transcript: an energy still dropping at `max_cycle` means "raise
`max_cycle` or use a second-order `algorithm`"; an *oscillating* energy means "add `damping`/`level_shift`,
or `smearing` for a small gap".
## Stability analysis
Convergence only means the SCF found *a* stationary point — not necessarily the **lowest** one.
`stability=True` tests whether the solution is a true minimum by checking the orbital-rotation Hessian for
negative eigenvalues (**internal** = a lower solution of the *same* reference; **external** = a lower solution
of a *broader* reference, e.g. RHF→UHF):
```python
h2 = qc.chk.new(atom="H 0 0 0; H 0 0 1.8", ao="cc-pvdz", unit="angstrom") # stretched H2
res = h2.scf(ref="r", stability=True).run()
dict(res.scf.stability)
# {'internal_stable': True, 'internal_eigenvalue': 0.3995,
# 'external_stable': False, 'external_eigenvalue': -0.1877, 'stable': False}
```
Stretched H₂ RHF is **externally unstable** (a negative external eigenvalue) — the closed-shell restriction
is wrong at long bond length, and a spin-broken **UHF** solution lies lower. This is the standard diagnostic
for bond dissociation and diradicals; the fix is to rerun as UHF, often with `guess(..., spin_break="mix")`
to break the symmetry.
## Worked example & exercise
```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
# HF vs a hybrid functional, same molecule and basis
for xc in (None, "b3lyp"):
s = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc=xc).run().scf
label = "HF" if xc is None else xc
print(f"{label:6} E = {s.energy:.6f} converged={s.converged} cycles={s.ncycle}")
# HF E = -76.026794 converged=True cycles=9
# b3lyp E = -76.420349 converged=True cycles=8
```
:::{exercise}
:label: ex-scf
1. You run the methyl radical CH₃ with `ref="r"`, `spin=2`. Which solver does qc-rs use, and is the result
spin-pure? What if you used `ref="u"`?
2. A UHF calculation gives $\langle S^2\rangle = 1.30$ for a doublet. Is that acceptable? What does it tell
you, and what would you try?
3. An RHF SCF on a stretched bond *converges* (`converged=True`) but you suspect it is not the ground state.
What one keyword checks this, and what result would confirm your suspicion?
:::
:::{solution} ex-scf
:class: dropdown
1. `ref="r"` on an **open** shell selects **ROHF** (see the reference table), which is **spin-pure**
($\langle S^2\rangle$ exact). `ref="u"` would give **UHF** — a lower energy but slightly spin-contaminated.
2. $\langle S^2\rangle = 1.30$ vs the exact doublet value $0.75$ is **heavy contamination** — the UHF
determinant is badly mixed with higher-spin states, so the energy and properties are suspect. Try **ROHF**
(spin-pure), or check whether the state is genuinely multireference (needs a correlated/multireference
method).
3. **`stability=True`**. A negative **external** eigenvalue (`external_stable: False`) confirms a lower
solution of a broader reference exists — rerun as UHF (with `spin_break="mix"`) to reach it.
:::
The SCF gives you a reference wavefunction and its energy. To recover the **correlation** it misses, the
[next chapter](post-scf.md) adds post-SCF methods — the RI-MP2 family.