# 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](../10-foundations/hartree-fock.md). 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](../00-intro/concepts.md) 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: $$ 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 . $$ $\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): ```python 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: ```python 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](../40-reference/checkpoints.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](scf.md). - **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](scf.md). ## 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: ```python 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} :label: ex-guess 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? ::: :::{solution} ex-guess :class: dropdown 1. **The same energy.** The SCF converges to the same fixed point regardless of the guess; `core` will just take more cycles to get there (it ignores electron–electron interaction entirely, so it starts further from the truth). 2. Converge **once** in a cheaper basis (say `cc-pvdz`) or at a nearby geometry, `save()` it, and pass it as `guess("read", source=...)` for the `cc-pvqz` runs — each expensive SCF then starts from an already near-converged density and needs far fewer of its costly cycles. 3. It **degrades `harris` to `gwh`**. `harris` (like `sad`/`sap`) is an all-electron atomic-superposition guess, which would assign the wrong valence electron count on an ECP atom; `gwh` builds from the valence-only ECP core Hamiltonian and is correct. ::: 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](scf.md).