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

Convergence: steering the SCF#

The SCF loop is iterative, so it needs an acceleration strategy to converge in a handful of cycles instead of hundreds. 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", runs DIIS (direct inversion of the iterative subspace — it extrapolates the next Fock matrix from previous cycles) 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`

The second-order strategies (qc/trah) build the exact orbital Hessian and take Newton steps — more work per cycle but far fewer cycles on a hard case:

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.