SCF: Hartree–Fock & KS-DFT

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: 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 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#

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). The three references differ in how they constrain the spin-orbitals (Part II):

`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:

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):
```
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).

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.

Algorithm 1 (CDIIS (Pulay) step)

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.

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).
Build the $n\times n$ Gram matrix $B_{ij} = \operatorname{Tr}\!\big(e_i^{\mathsf T} e_j\big)$.
Solve the Lagrange-bordered linear system for the extrapolation coefficients $\mathbf c$ (constrained to $\sum_i c_i = 1$, enforced by the multiplier $\lambda$): $$

(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}. $$
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 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 $$

(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, $$ 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.

        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<br/>energy-model step"]
    Q -->|"small (near)"| CD["CDIIS<br/>Pulay extrapolation"]
    Q -->|"stalling"| SO["2nd-order finisher<br/>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.

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, 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:

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.

Reading and diagnosing a run#

After .run(), the scf accessor exposes the result; run(log=...) (from the quickstart) streams the live cycle-by-cycle transcript that lets you watch convergence:

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):

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#

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 4

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"?
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?
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 to Exercise 4

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.
$\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).
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 adds post-SCF methods — the RI-MP2 family.