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.

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.