The initial guess

The initial guess#

You now have a molecule and a basis. Before the SCF can run, it needs somewhere to start — an approximate density to build the very first Fock matrix from. That starting point is the initial guess. Most of the time you never think about it: qc-rs picks a good default automatically. But understanding it pays off exactly when a calculation is hard to converge — so this short chapter explains what the guess is, why it (usually) does not change your answer, and how to change it when you need to.

Theory: why an SCF needs a guess#

Recall the self-consistent-field loop from Hartree–Fock. The Fock matrix \(\mathbf F\) is built from the density \(\mathbf D\) (through the Coulomb and exchange operators \(\mathbf J[\mathbf D]\), \(\mathbf K[\mathbf D]\)), but the density comes from diagonalizing \(\mathbf F\). That circularity —

\[ \mathbf D \;\longrightarrow\; \mathbf F[\mathbf D] \;\xrightarrow{\;\mathbf{FC}=\mathbf{SC}\boldsymbol\varepsilon\;}\; \mathbf C \;\longrightarrow\; \mathbf D' \;\longrightarrow\; \cdots \]

is why the SCF is iterative. But it also means the loop cannot take its first step without a density already in hand. The initial guess provides that starting \(\mathbf D_0\) (or a starting set of orbitals).

Here is the key idea. The SCF converges to a fixed point — the self-consistent density — and for a well-behaved system that fixed point does not depend on where you started. So:

Important

The guess changes the path, not the destination A different initial guess gives the same converged energy; it only changes how many cycles the SCF needs (and, for difficult systems, whether it converges at all). A guess close to the final density converges fast; a crude one takes more cycles. You do not trade accuracy for a cheaper guess — only speed and robustness.

We will see this directly in the worked example below: six different guesses for water all land on \(-76.026794\,E_h\), taking between 8 and 12 cycles.

The default: SAD#

qc-rs’s default is the superposition of atomic densities (SAD). The idea is physical and cheap: a molecule’s density is, to first approximation, just its atoms’ densities overlaid. SAD precomputes a spherically-averaged density for each element (from a small atomic SCF) and sums them at the nuclear positions to form \(\mathbf D_0\). It needs no molecular integrals beyond what it already has, is robust across the periodic table, and is usually close enough to converge in a handful of cycles.

You do not have to ask for it. As Core concepts described, run() auto-inserts a sad guess when an SCF needs a starting state and none exists. (If an atomic density cannot be built — an untabulated element, or an ECP-only heavy atom — SAD falls back to gwh, then to the bare core Hamiltonian.)

The menu of guesses#

You can override the default with guess("..."). The available guesses, from most physical to crudest, with the cycle count each took for water/cc-pvdz (RHF) in the example below:

guess(...)

what it does

cycles*

"sad" (default)

superposition of atomic densities — robust, physical

9

"minao"

fixed minimal-AO occupation density (PySCF minao); no atomic SCF, so cheaper than SAD

9

"harris"

Harris-functional density (a non-self-consistent superposition)

8

"sap"

superposition of atomic potentials (\(\mathbf F = \mathbf H_{\text{core}} + \mathbf J[\mathbf D_{\text{SAD}}]\), Coulomb only)

10

"gwh"

generalized Wolfsberg–Helmholz (from the core Hamiltonian); ECP-aware

12

"core"

bare core Hamiltonian \(\mathbf T + \mathbf V\) — ignores all electron–electron interaction; crudest

11

"read"

import/project orbitals from another checkpoint (restart)

—

*Cycle counts are system-dependent — do not read this as a universal ranking. For water they are close; for a hard system the spread (and which guess wins) can be very different. The point is that they all reach the same energy.

Theory: what gwh and sap actually compute#

core and gwh are the two guesses that need no two-electron integrals at all — useful when even building \(\mathbf J,\mathbf K\) once is unwelcome, or as the universal ECP-safe fallback (see below). core simply diagonalizes the bare one-electron Hamiltonian \(\mathbf H_{\text{core}}=\mathbf T+\mathbf V\) — correct in form, but it ignores electron–electron repulsion entirely, so its orbitals are a poor starting point.

GWH (generalized Wolfsberg–Helmholtz) patches this with a semi-empirical off-diagonal correction built purely from the overlap matrix and a table of atomic orbital energies — still no two-electron integrals:

\[\begin{split} H_{\mu\nu}^{\mathrm{GWH}} = \begin{cases} \epsilon_\mu & \mu=\nu, \\[4pt] K_{\mathrm{GWH}}\, S_{\mu\nu}\, \dfrac{\epsilon_\mu+\epsilon_\nu}{2} & \mu\ne\nu, \end{cases} \qquad \mathbf H^{\mathrm{GWH}}\mathbf C=\mathbf S\mathbf C\boldsymbol\varepsilon . \end{split}\]

\(\epsilon_\mu\) is a per-AO valence orbital energy from a tabulated element-by-angular-momentum table (H–Ar); an AO the table does not cover falls back to the bare diagonal \(H^{\text{core}}_{\mu\mu}\). qc-rs fixes \(K_{\mathrm{GWH}}=1.75\) (the literature range is \(1.5\)–\(2.0\)). The intuition: two AOs that overlap strongly and sit at similar energy should couple strongly — exactly what this off-diagonal term builds — which gives the eigenvectors a first hint of bonding the bare core Hamiltonian cannot see.

SAP (superposition of atomic potentials) goes one step further and folds in the electron repulsion from the SAD density — Coulomb only, no exchange:

\[ \mathbf F^{\mathrm{SAP}} = \mathbf H_{\text{core}} + \mathbf J[\mathbf D_{\text{SAD}}]. \]

Because \(\mathbf J[\cdot]\) is linear in the density, \(\mathbf J[\mathbf D_{\text{SAD}}] = \sum_A \mathbf J[\mathbf D_A]\) — the Fock-like operator is exactly a sum of each atom’s own Hartree potential, which is what “superposition of atomic potentials” names. It costs one Coulomb build (cheaper than a full SCF cycle, since no exchange and no self-consistency), buying a guess close enough to shave a cycle or two off DIIS (10 vs. 12 for gwh in the table above).

Usage#

The common case is doing nothing — let run() insert SAD. To override, add a guess step before the SCF (functional or method-chain form, like every other step):

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

# explicit guess (method-chain)
mychk = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").guess("gwh").scf(ref="r").run()

# equivalently, functional form
base  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom")
mychk = qc.scf(qc.guess(base, "gwh"), ref="r").run()

Restarting from a previous result: guess("read")#

guess("read", source=...) imports the orbitals from another checkpoint, projecting them into the current basis. Two everyday uses:

  • Restart a calculation from a saved .qch5 — continue, or re-run a property on an existing SCF.

  • Basis stepping — converge cheaply in a small basis, then project that result as the guess for the big-basis SCF, which then needs far fewer expensive cycles:

small = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
small.save("water_dz.qch5")

big = qc.chk.new(atom=water, ao="cc-pvtz", unit="angstrom") \
        .guess("read", source="water_dz.qch5").scf(ref="r").run()

The read details (MO projection, symmetry/irrep handling) are in the checkpoint reference, qc-chk.md.

When the guess actually matters#

For routine closed-shell molecules the default just works. The guess earns your attention in three cases:

  • Hard-to-converge systems — transition metals, near-degeneracies, diradicals. A better guess (or a read from a related calculation) can be the difference between smooth convergence and an SCF that oscillates. Pair it with the convergence tools in the SCF chapter.

  • ECP basis sets — the atomic-superposition guesses (sad/sap/harris) are all-electron and would put the wrong number of valence electrons on an ECP atom, so qc-rs automatically degrades them to gwh (which diagonalizes the valence core Hamiltonian \(\mathbf T + \mathbf V_{\text{nuc}}[Z_{\text{eff}}] + \mathbf V^{\text{ECP}}\)). gwh and core are ECP-correct directly.

  • Broken symmetry & excited states — a UHF diradical often needs guess(..., spin_break="mix") to break spatial symmetry, and a ΔSCF excited state needs a non-Aufbau occupation (guess(..., occupation=...) with mom=True). These are covered in the SCF chapter.

Worked example: same energy, different number of cycles#

This is the whole idea of the chapter in one script — every guess converges to the identical energy, but the count of SCF cycles varies:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

for g in ("sad", "minao", "harris", "sap", "core", "gwh"):
    chk = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").guess(g).scf(ref="r").run()
    print(f"{g:8}  E = {chk.scf.energy:.6f}   cycles = {chk.scf.ncycle}")

# sad       E = -76.026794   cycles = 9
# minao     E = -76.026794   cycles = 9
# harris    E = -76.026794   cycles = 8
# sap       E = -76.026794   cycles = 10
# core      E = -76.026794   cycles = 11
# gwh       E = -76.026794   cycles = 12

The energies are bit-for-bit identical; only the cycles column moves. That is the guess doing its one job: setting the starting point, not the answer.

Exercise 3

  1. Without running it, what converged RHF energy do you expect from guess("core") versus guess("sad") for the same molecule and basis? Why?

  2. You must run 50 SCFs on the same molecule at cc-pvqz (an expensive basis). Sketch a strategy using guess("read") to make each one converge in fewer cycles.

  3. You give a heavy-metal complex an ECP basis and request guess("harris"). What does qc-rs actually use, and why?

With the molecule, basis, and starting point settled, we can finally turn to the SCF itself — the references, functionals, and the convergence toolkit — in the next chapter.