Initial guesses¶

Every SCF needs a starting point. The initial guess supplies the first density matrix (or, equivalently, the first set of orbital coefficients) that the SCF iteration refines. A good guess can cut iteration counts in half and steer SCF away from spurious local minima; a bad one can stall the iteration entirely or, in periodic calculations, drive the energy off to nonsense values (O(10⁴) Ha) before it ever recovers.

vibe-qc uses one unified engine for initial-guess construction, vibeqc::GuessEngine, shared between every molecular SCF driver (RHF / UHF / RKS / UKS) and every Gamma periodic SCF driver (C++ + the Python Ewald / GDF families). InitialGuess selects the method; the AUTO default picks one for you based on the system.

Quick start¶

import vibeqc as vq
from vibeqc import Atom, Molecule, BasisSet, RHFOptions, run_rhf

mol = Molecule([
    Atom(8, [0.0,  0.0,  0.0]),
    Atom(1, [0.0,  1.43, -0.98]),
    Atom(1, [0.0, -1.43, -0.98]),
])
basis = BasisSet(mol, "6-31g*")

opts = RHFOptions()
opts.initial_guess = vq.InitialGuess.AUTO     # (default) picks SAP for this case

result = run_rhf(mol, basis, opts)

The nine available kinds:

vq.InitialGuess.AUTO       # default: engine picks per system
vq.InitialGuess.HCORE      # diagonalise the core Hamiltonian
vq.InitialGuess.SAD        # superposition of atomic densities
vq.InitialGuess.SAP        # superposition of atomic potentials (Lehtola 2020)
vq.InitialGuess.PATOM      # SAD + one in-field re-polarisation step
vq.InitialGuess.HUECKEL    # parameter-free generalised Wolfsberg-Helmholz
vq.InitialGuess.MINAO      # minimal-AO (ANO-RCC) reference projection
vq.InitialGuess.READ       # restart from a prior result / .qvf / Molden file
vq.InitialGuess.FRAGMO     # superposition of converged fragment densities

All nine are implemented today. SAP, PATOM, HUECKEL, and MINAO additionally work for periodic closed-shell SCF (Gamma and multi-k RHF/RKS; see their sections). SAP, PATOM, HUECKEL, and MINAO also reach the Python periodic UHF/UKS drivers. READ works for molecular, Gamma periodic, and closed-shell multi-k GDF/GPW/GAPW restarts from in-memory results or current .qvf archives (see the READ section). FRAGMO is molecular-only: a periodic SCF selecting it raises a clear error because a crystal has no finite-fragment partition.

What an initial guess actually does¶

The SCF iteration solves

\[ \mathbf{F}(\mathbf{D})\,\mathbf{C} = \mathbf{S}\,\mathbf{C}\,\boldsymbol{\varepsilon}, \qquad \mathbf{D} = 2\,\mathbf{C}_{\rm occ}\,\mathbf{C}_{\rm occ}^{\!T}, \]

a fixed-point problem where the Fock matrix $\mathbf{F}$ depends on the density $\mathbf{D}$, which is built from the occupied orbitals $\mathbf{C}_{\rm occ}$, which come from diagonalising $\mathbf{F}$. To start the loop you need either an initial $\mathbf{D}$ (a density-mode guess), or an initial $\mathbf{F}_{\rm guess}$ that the driver can diagonalise once to get $\mathbf{C}$ and then $\mathbf{D}$ (a Fock-mode guess), or an initial $\mathbf{C}$ directly (an MO-mode guess from a prior calculation).

The literature classifies guess methods by which of those three artifacts they produce:

Mode	Methods	Engine output field
Density	`SAD`, `PATOM`, `MINAO`	`D`
Fock	`HCORE`, `SAP`, `HUECKEL`	`F_guess`

READ and FRAGMO sit outside this table: both are resolved in the run_* wrapper rather than the engine, injecting a density as opts.read_density (closed shell) or opts.read_density_alpha / opts.read_density_beta (open shell). READ injects a prior density; FRAGMO injects a freshly assembled block-diagonal superposition of converged fragment densities. The READ and FRAGMO sections below cover them.

The GuessEngine returns a tagged result and each SCF driver consumes whichever artifact is populated. Closed-shell drivers see GuessClosedShellResult; open-shell drivers see GuessOpenShellResult, which carries per-spin matrices.

Initial guesses cannot change the converged minimum. Two different guesses on a well-behaved SCF reach the same total energy to numerical precision (typically below $10^{-10}$ Ha). What they change is the path: iteration count, intermediate energies, and, for systems with multiple stationary points, which stationary point the SCF lands on. The OH/6-31G* UHF false-minimum trap is the canonical example: HCore converges to a local min ~0.16 Ha above the true UHF minimum; SAD reaches the right minimum.

Multiplicity and the initial guess¶

Multiplicity and the initial guess do different jobs. The multiplicity fixes how many electrons of each spin the SCF must converge to; the guess supplies a starting density that respects those counts but does not, on its own, decide where the spin sits.

Multiplicity sets the spin counts. From the multiplicity $2S+1$ and the total electron count $n_e$ (set by the nuclear charges and the molecular or cell charge), each driver derives the alpha and beta occupations with the usual convention

\[ n_\alpha = \frac{n_e + (\text{mult}-1)}{2}, \qquad n_\beta = \frac{n_e - (\text{mult}-1)}{2}, \]

so $n_\alpha - n_\beta = \text{mult}-1 = 2S$ is the number of net unpaired electrons. That is the only spin quantity you set directly. vibe-qc’s open-shell SCF is unrestricted (UHF / UKS) with these two counts held fixed; $\langle S^2\rangle$ is not constrained, so the converged value is whatever the unrestricted solution yields (hence the spin-contamination check in the worked examples).

The guess builds to those counts; its spin symmetry follows a published, deliberate policy. Given $n_\alpha$ and $n_\beta$, the open-shell guess produces a starting point with the right counts and a spin symmetry decided by the electron counts and the requested seed, not by atomic Hund rules alone:

$n_\alpha = n_\beta$ (spin-balanced): spin-averaged atomic occupations – each atom contributes $\tfrac12\,f_i$ to both spins from the spherically-averaged atomic SCF, so $\mathbf{D}_\alpha = \mathbf{D}_\beta$ exactly. This is the deterministic symmetric start (the Van Lenthe 2006 standard: the UHF start comes from the closed-shell-type density), and it keeps a genuine closed-shell singlet bit-symmetric through the whole SCF.
$n_\alpha \ne n_\beta$: Hund-split atomic densities – each atom contributes its natural majority-alpha configuration ($f_i^\alpha$, $f_i^\beta$ from spin_resolve_occupations), so the starting spin density is per-atom polarised. This is what lands FeCl$_3$-class open shells in the right basin (BUG 88): a spherical spin density there converges to an excited solution.
The UKS packaging diagonalises $\mathbf{F}_{\rm SAD}$ once from the total density and keeps the first $n_\alpha$ and $n_\beta$ columns.

The physical spin density develops during SCF. For a single radical (one unpaired electron) there is only one way to place the spin and the SCF finds it. For a system with several competing spin solutions the default symmetric start biases toward the most symmetric one; the default-on internal stability analysis (Seeger-Pople 1977; Lehtola 2020 section 10) then detects a saddle and attempts to descend toward a broken-symmetry minimum – verified on H2 past the Coulson-Fischer point, ozone, twisted ethylene, p-benzyne, and Cr2. The escape searches both signs of the unstable orbital rotation and rechecks the final Hessian before accepting the result. The first four systems reach identical solutions before and after the spin-averaging fix (handovers/HANDOVER_GUESS_SPIN_SYMMETRY.md). A deliberate antiferromagnetic / ferrimagnetic pattern is requested explicitly via atomic_spins (below), never imposed by the default guess.

Seeding a broken-symmetry spin pattern (the CRYSTAL ATOMSPIN analogue). Molecular UHF and UKS accept atomic_spins, a per-atom list of +1 / -1 / 0 in atom order (+1 = majority alpha, -1 = majority beta, 0 = unpolarised). It seeds a broken-symmetry SAD start in which each tagged atom contributes a Hund’s-rule spin-resolved atomic density, assembled block-diagonally into $\mathbf{D}_\alpha$ / $\mathbf{D}_\beta$:

opts = UHFOptions()
opts.atomic_spins = [+1, -1]   # atom 0 spin-up, atom 1 spin-down (antiparallel)

This is what lets you initialise an antiferromagnetic or ferrimagnetic state that the spin-symmetric SAD split cannot reach (think alternating metal-site spins in an oxide). The SCF then refines the pattern while the total $n_\alpha - n_\beta$ stays fixed by the multiplicity. atomic_spins is consulted only when the guess resolves to SAD (including AUTO on an open shell); pairing it with another guess, or giving a wrong-length list, raises.

atomic_spins works for periodic UHF/UKS too: the per-atom tags seed the g=0 cell block (in unit-cell atom order), which Bloch-sums to a broken-symmetry $\mathbf{D}(\mathbf{k})$, the route to an AFM/ferrimagnetic periodic start. It is wired on the Γ UHF/UKS Ewald, GDF, and BIPOLE drivers and the multi-k UKS Ewald driver (the flagship for AFM solids); the multi-k UHF Ewald driver, which has no broken-symmetry guess hook, raises. Pass it through vq.run_periodic_job(..., method="UHF"/"UKS", atomic_spins=[...]), or set opts.atomic_spins on a driver directly.

Holding the pattern through early iterations (SPINLOCK). A broken-symmetry start can still wash out to the symmetric solution in the first few iterations on a hard system. SPINLOCK (the CRYSTAL analogue) holds it, in either of two modes on molecular UHF / UKS, set through opts.spinlock_mode, opts.spinlock_iterations, and (for the schedule) opts.spinlock_value:

SpinlockMode.PATTERN_HOLD keeps the seeded broken-symmetry occupied set in place by maximum-overlap selection (MOM, Gilbert/Besley/Gill 2008) instead of plain Aufbau for the first spinlock_iterations cycles, then releases. Pair it with atomic_spins to stop an antiferromagnetic seed from collapsing early.
SpinlockMode.SPIN_SCHEDULE runs a two-phase SCF: it converges at a locked $n_\alpha - n_\beta = $ spinlock_value for the first spinlock_iterations cycles, then restarts from that density at the multiplicity target. This is CRYSTAL’s SPINLOCK n nstep; use it to drive the SCF through a high-spin intermediate into the state you want.

opts = UHFOptions()
opts.atomic_spins = [+1, -1]                      # seed the broken-symmetry start
opts.spinlock_mode = vq.SpinlockMode.PATTERN_HOLD
opts.spinlock_iterations = 8                      # hold 8 cycles, then release

Periodic SPINLOCK. Both modes also work for periodic open-shell SCF: PATTERN_HOLD on the Γ UHF/UKS Ewald, GDF and BIPOLE drivers and the multi-k UKS Ewald driver (per-k per-spin MOM hold, the path that protects an atomic_spins AFM seed against collapse on multi-k); SPIN_SCHEDULE on the Γ UHF/UKS Ewald, Γ GDF and Γ BIPOLE drivers (a two-phase SCF that locks a high-spin state, then restarts at the target multiplicity from that density via the READ machinery). Both are reachable from the high-level runner: vq.run_periodic_job(..., method="UHF"/"UKS", spinlock="pattern_hold" / "spin_schedule", spinlock_iterations=N[, spinlock_value=M]). Periodic SPINLOCK on a driver that does not implement the requested mode (the multi-k UHF Ewald path for any mode, or SPIN_SCHEDULE on the multi-k UKS Ewald and BIPOLE paths) fails closed with a clear error rather than silently ignoring the request.

On every driver that implements it (molecular UHF/UKS, Γ UHF/UKS Ewald, GDF, BIPOLE, and multi-k UKS Ewald), PATTERN_HOLD also suspends the SCF accelerator (plain DIIS as well as the EDIIS / ADIIS / KDIIS variants) for the hold window (density damping stays live) and starts the extrapolation history fresh at release. MOM only selects which orbitals are occupied; it cannot stop Fock extrapolation from rotating the occupied orbitals themselves toward the spin-symmetric solution, so letting the accelerator extrapolate across held iterations would quietly collapse the seed inside the hold window: exactly the failure the hold exists to prevent.

Still pending: BETALOCK (locking only the beta count), SPIN_SCHEDULE on the multi-k routes, and per-element starting occupations; those are on the roadmap.

Two further occupation-related controls ship, though neither seeds a spin pattern the way atomic_spins does:

Control	Where	What it does	What it is not
`use_mom=True` (maximum overlap method)	periodic multi-k RHF / RKS	Holds the occupied set fixed by tracking maximum orbital overlap across iterations instead of refilling by energy, so a chosen non-Aufbau occupation survives the SCF.	A way to specify a per-site spin pattern from scratch; it follows the orbitals you start from.
`smearing_temperature = T` (Fermi-Dirac)	periodic RHF / RKS	Replaces hard Aufbau with fractional occupations near $E_F$, easing convergence of metals and near-degenerate states.	A lock; it smooths occupations rather than fixing a pattern.

The underlying methods are MOM (Gilbert, Besley, Gill 2008) and the Mermin finite-temperature functional (Mermin 1965); both are in the references.

The methods, with math¶

HCORE, core-Hamiltonian guess¶

The cheapest possible guess. Set the initial Fock to the one-electron core Hamiltonian and ignore electron-electron repulsion entirely:

\[ \mathbf{F}_{\rm HCORE} \;=\; \mathbf{T} + \mathbf{V}_{\rm ne} \]

where $\mathbf{T}$ is the kinetic-energy matrix and $\mathbf{V}_{\rm ne}$ is the nuclear-attraction matrix. Diagonalise once, build $\mathbf{D}$ from the occupied orbitals, hand off to the SCF loop.

Cost: trivial, no extra integrals beyond those the driver needs anyway.

Pathology: the shell structure is wrong. Without screening from the inner electrons, all the orbitals localise on the heaviest atom and the outer-valence orbitals have the wrong sign of energy. On molecular systems with light atoms this is a mild slowdown (2-3 extra iterations); on ionic insulators with deep cores the first Fock build from $\mathbf{D}_{\rm HCORE}$ can push the energy into the $\pm 10^4$ Ha range, which DIIS then has to climb back out of, the NaCl-bombing failure mode.

Use HCORE explicitly when you need a deterministic, system- independent reference (e.g. for unit tests that compare implementations).

SAD, superposition of atomic densities¶

Build the molecular density as a block-diagonal sum of converged atomic densities. For each unique element $A$ in the molecule, run a small spherically-averaged atomic SCF with fractional occupations:

\[ \mathbf{D}_{A,\rm atom} \;=\; \sum_{i} f_i \, \mathbf{C}_i^{(A)} \mathbf{C}_i^{(A)\,T} \]

where $f_i \in [0, 2]$ are the orbital occupations. Because the spherically-averaged atomic Fock is block-diagonal in the angular momentum $\ell$, every atomic MO is a pure-$\ell$ function, and the occupations are filled per $\ell$-channel up to the element’s ground-state electron count for that channel (Aufbau within each channel, with degenerate orbitals sharing equal fractional occupation so the result stays spherical). The per-$\ell$ counts come from the neutral-atom ground-state configuration, the Madelung order plus the standard anomalies (Cr, Cu, Pd, the lanthanide/actinide irregulars, …), the same data PySCF, Psi4, and ORCA pin for their atomic guesses. This matters for the transition metals: a single global eigenvalue sort over all orbitals mis-occupies the 3d row, because the 3d/4s/4p near-degeneracy makes the bare-atom orbital order disagree with the physical configuration (it would seat Ni at 3d¹⁰ 4s⁰ rather than the physical 3d⁸ 4s²). Pinning per channel reproduces the experimental d-shell filling for all of Sc-Zn. Then assemble the molecular density by placing each atomic block on the diagonal:

\[ \mathbf{D}_{\rm SAD} \;=\; \bigoplus_{A \in \text{atoms}} \mathbf{D}_{A,\rm atom}. \]

For periodic systems the engine places $\mathbf{D}_{\rm SAD}$ at the $g=0$ cell and leaves other cells zero; the SCF builds the rest by Bloch summation on the first iteration.

Per-spin packaging. Open-shell drivers see two flavours:

UHF (proportional split): $\mathbf{D}_{\alpha} = (n_\alpha / n_e)\,\mathbf{D}_{\rm SAD}$, $\mathbf{D}_{\beta} = (n_\beta / n_e)\,\mathbf{D}_{\rm SAD}$. Cheap; spin polarisation develops on the first iteration via $\mathbf{F}_\alpha \ne \mathbf{F}_\beta$ when $n_\alpha \ne n_\beta$.
UKS (Fock-diagonalise split): build the closed-shell-style $\mathbf{F}_{\rm SAD} = \mathbf{H}_{\rm core} + \mathbf{J}(\mathbf{D}_{\rm SAD}) - \tfrac{1}{2}\,\mathbf{K}(\mathbf{D}_{\rm SAD})$, diagonalise it, and split the occupied columns into per-spin densities. Costs one extra Fock build at SCF start but lands closer to the Kohn-Sham minimum.

The engine selects the packaging automatically based on whether a JKBuilder* is supplied (UKS supplies one; UHF does not).

Cost: one atomic SCF per unique element (cached). For typical molecules this is microseconds; even on a heavy atom (e.g. Z=86 Rn) the atomic SCF is fast because it’s spherically symmetric and a single S-shell expansion.

Strength: robust on every system class vibe-qc supports. SAD is the right default for ionic insulators, transition-metal complexes, and any system where the Hcore failure mode would bite.

SAP, superposition of atomic potentials (Lehtola 2020)¶

A Fock-mode method designed by Lehtola, Visscher, and Engel to give shell-structure-correct starting orbitals without SAD’s atomic- SCF cost. Each atom $A$ contributes an erf-screened Coulomb potential fit to a numerically-exact atomic SCF:

\[ V_{\rm SAP}^{(A)}(r) \;=\; -\sum_{i} c_i^{(A)}\,\frac{\operatorname{erf}\!\left(\sqrt{\alpha_i^{(A)}}\,r\right)}{r} \]

where $(\alpha_i^{(A)}, c_i^{(A)})$ are tabulated for $Z = 1$ to $86$ (see bundled datasets below). The total SAP potential is

\[ V_{\rm SAP}(\mathbf{r}) \;=\; \sum_{A} V_{\rm SAP}^{(A)}(|\mathbf{r} - \mathbf{R}_A|) \]

and the initial Fock matrix is

\[ \mathbf{F}_{\rm SAP} \;=\; \mathbf{T} + \mathbf{V}_{\rm SAP}. \]

Note that $\mathbf{V}_{\rm SAP}$ replaces $\mathbf{V}_{\rm ne}$ entirely, the screened potential already includes the bare nuclear attraction in its $r \to 0$ limit, and adding $\mathbf{V}_{\rm ne}$ would double-count it.

The matrix element of one $(\alpha_i, R_A)$ primitive is exactly libint2’s Operator::erf_nuclear with attenuation $\omega = \sqrt{\alpha_i}$ and a “point charge” $-c_i$ at $R_A$:

\[ \langle\chi_\mu \mid V_{\rm Gaussian}(\alpha, R_A) \mid \chi_\nu\rangle \;=\; \langle\chi_\mu \mid \frac{\operatorname{erf}(\sqrt{\alpha}\,|\mathbf{r}-\mathbf{R}_A|)}{|\mathbf{r}-\mathbf{R}_A|} \mid \chi_\nu\rangle \]

So the engine loops over (atom $A$, primitive $i$) pairs, calls libint once per pair, and accumulates into $\mathbf{V}_{\rm SAP}$.

Cost: one libint erf_nuclear engine call per (atom, primitive) pair, typically 5-30 primitives per atom × $n_{\rm atoms}$, then one diagonalisation. No atomic SCF; no Fock build. Cheaper than SAD on heavy elements where the atomic SCF dominates.

Strength: unlike HCore, the shell structure is correct, the first iteration sees orbitals that look like Hartree-Fock orbitals rather than artefacts of nuclear attraction alone. Iteration counts on light-atom molecules typically match SAD within ±1; on heavy atoms SAP wins by 30-50% per [Lehtola 2020] Tables I-II.

Limitations:

Periodic SAP is wired for closed-shell RHF/RKS and Python UHF/UKS. The lattice-summed SAP potential (compute_vsap_lattice, built by an Ewald split, validated against the molecular SAP matrix) gives the Fock-mode guess F_SAP = T + V_SAP, diagonalised at Γ or each k. Supported by the Γ RHF driver (run_rhf_periodic_gamma) and the multi-k RHF/RKS drivers (run_rhf_periodic / run_rks_periodic), plus the Python Gamma/multi-k Ewald, GDF, and BIPOLE UHF/UKS drivers. The open-shell Python routes Gamma-fold the same SAP Fock and fill the requested alpha/beta occupations; AUTO still resolves to SAD for periodic.
Element coverage. The bundled tables cover Z = 1-86. Elements outside the table fall back to bare nuclear attraction (-Z/r), so SAP gracefully degrades to “SAP for known atoms, Hcore for the rest” on exotic systems.

SAP datasets¶

Two atomic-potential libraries ship with vibe-qc in python/vibeqc/basis_library/basis/:

File	Reference	Default for
`sap_helfem_large.g94`	All-electron, Helfem fully-numerical atomic SCF	Non-relativistic calculations
`sap_grasp_large.g94`	Relativistic Dirac-Hartree-Fock, GRASP	x2c / DKH calculations

The engine selects sap_helfem_large by default; the GRASP table becomes the AUTO choice once relativistic-Hamiltonian detection lands (the data is already in the wheel; only the selection heuristic is missing).

Where the tables are looked up. SAP tables resolve at run time, with the same precedence libint uses for orbital basis sets: $LIBINT_DATA_PATH/basis/<name>.g94 first, then the data directory recorded when the binary was built. import vibeqc already points $LIBINT_DATA_PATH at the bundled library, so a wheel installed anywhere – or a relocatable bundle whose build tree is long gone – finds its tables without you setting anything. Set LIBINT_DATA_PATH yourself only to point SAP (and the orbital bases) at a different library. If a table is genuinely missing, the error lists every path that was searched.

Citation requirement. Any published result that uses the SAP initial guess must cite

S. Lehtola, L. Visscher, E. Engel, Efficient Implementation of the Superposition of Atomic Potentials Initial Guess for Electronic Structure Calculations in Gaussian Basis Sets, J. Chem. Phys. 152, 144105 (2020).

See docs/license.md §3a for the full bundled-data inventory.

AUTO, pick the right guess for me¶

The default for all six SCF options structs. AUTO inspects the molecule and chosen system flags and resolves to a concrete kind via GuessEngine::resolve_auto:

System	AUTO resolves to	Reason
Periodic, any system	SAD	AUTO picks SAD for robustness on ionic insulators (the NaCl-bombing fix); explicit `SAP`, `PATOM`, `HUECKEL`, and `MINAO` are wired for closed-shell RHF/RKS and Python UHF/UKS
Molecular, open-shell	SAD	Spin polarisation develops via per-spin Fock asymmetry
Molecular, transition or f-block atom	SAD	SAP closed-shell averaging can land the wrong d-shell occupation; select `PATOM` (now shipped, molecular) for a re-polarised SAD start
Molecular, closed-shell, light atoms	SAP	Shell-structure-correct, no atomic SCF cost

The hints (is_open_shell, has_transition_metal, is_periodic) are filled in automatically by each SCF wrapper:

is_open_shell from n_alpha != n_beta
has_transition_metal from a Z-range scan over the molecule (Sc-Zn, Y-Cd, La/Hf-Hg, Ce-Lu, Ac/Rf-Cn, Th-Lr)
is_periodic from which entry point the caller used

The resolved kind is recorded on the result so the choice is visible in the SCF log:

result = run_rhf(mol, basis, opts)
# (look in the SCF log; "initial guess: SAP (Lehtola/Visscher/Engel 2020 …)"

PATOM, HUECKEL, and MINAO¶

Three further guesses landed together. PATOM and molecular MINAO are density-mode guesses that return a starting density the SCF consumes directly; HUECKEL is Fock-mode and is diagonalised once to build the starting density. For periodic closed-shell SCF, PATOM, HUECKEL, and MINAO have lattice-specific Gamma and multi-k RHF/RKS implementations. Python periodic drivers also pass the periodic context needed by SAP, the HUECKEL lattice Fock helper, the MINAO density projection, and PATOM’s route-local in-field Fock step. Open-shell Python SAP/HUECKEL build per-spin densities from those helpers; open-shell MINAO uses the projected total density with the proportional spin split; open-shell PATOM starts from SAD and applies one HF-like periodic in-field re-polarisation step per spin.

PATOM takes the SAD density and applies one in-field mean-field re-polarisation step, letting the superposed free-atom densities relax in the molecular field before the SCF proper. It is the density-mode upgrade of SAD for cases where SAD’s spherical atoms are too rigid, notably transition-metal d-shell ordering. run_uhf forwards its JKBuilder to the engine only for PATOM, so no other guess’s packaging changes. The periodic closed-shell drivers run the same one-step re-polarisation: Gamma forwards the periodic JKBuilder, while multi-k builds the in-field Fock in real space and diagonalises it at each k-point. The Python Gamma/multi-k Ewald, GDF, and BIPOLE UHF/UKS drivers use their own periodic J/K builders for the open-shell per-spin in-field step. (ORCA PAtom; builds on Van Lenthe SAD.)
HUECKEL is a parameter-free generalised Wolfsberg-Helmholz (extended-Hückel) guess. The off-diagonal elements are $H_{\mu\nu} = \tfrac{1}{2} K\, S_{\mu\nu}\,(H_{\mu\mu} + H_{\nu\nu})$ with $K = 1.75$, built over computed atomic-orbital energies rather than tabulated parameters, then diagonalised to a density. (Wolfsberg and Helmholz 1952; Hoffmann 1963; the parameter-free variant of Lehtola 2019.) For periodic closed-shell RHF/RKS the same GWH form is built as a lattice Fock, $H_{\mu\nu}(\mathbf{g})$, from the lattice overlap $S_{\mu\nu}(\mathbf{g})$ and diagonalised at Gamma or each k-point. The closed-shell Python periodic drivers Gamma-fold that lattice Fock and inject the resulting g=0 density. Accepts initial_guess="huckel" or "hueckel".
MINAO builds free-atom reference densities in the ANO-RCC minimal contraction and projects them onto the target basis (Knizia 2013 minimal-AO concept; ANO-RCC reference data). For periodic systems the projection is lattice-summed, P = S(Γ)^{-1} S_tm(Γ) with the cross-basis overlap S_tm(Γ) = Σ_g <χ(0)|χ_ref(g)>; it is wired into the closed-shell Γ and multi-k RHF/RKS drivers and the closed-shell Python periodic density seam, seeding the g=0 density cell (like SAD).

FRAGMO, superposition of converged fragment densities (molecular)¶

Where SAD superposes converged atomic densities, FRAGMO superposes converged fragment densities. You partition the molecule into fragments; each fragment is converged with its own SCF, and the fragments’ density matrices are assembled block-diagonally into the supersystem guess:

\[ \mathbf{D}_{\rm FRAGMO} \;=\; \bigoplus_{F \in \text{fragments}} \mathbf{D}_{F}, \]

each fragment block $\mathbf{D}_F$ placed in the supersystem AO rows/columns of that fragment’s atoms (intra-fragment cross terms preserved; inter-fragment blocks left at zero). FRAGMO is resolved in the run_* wrapper and injected via opts.read_density exactly like READ, so one path covers RHF / UHF / RKS / UKS.

import vibeqc as vq
from vibeqc import Fragment, run_rhf

opts = vq.RHFOptions()
opts.initial_guess = vq.InitialGuess.FRAGMO
# Two fragments: atoms 0-2 and atoms 3-5 (each neutral closed-shell).
result = run_rhf(mol, basis, opts, fragments=[[0, 1, 2], [3, 4, 5]])

# A charged / open-shell fragment uses an explicit Fragment spec:
result = run_rhf(mol, basis, opts,
                 fragments=[Fragment(atoms=[0, 1, 2], charge=+1),
                            Fragment(atoms=[3, 4, 5], charge=-1)])

Each fragment is converged with the supersystem’s method and functional (a closed-shell fragment with RHF/RKS, an odd-electron or multiplicity > 1 fragment with UHF/UKS); for an open-shell supersystem each fragment’s α/β blocks feed the matching supersystem spin.

Strength. For weakly-interacting assemblies (hydrogen-bonded dimers, van-der-Waals complexes, host-guest systems) the fragments already carry almost all of the electronic structure, so the guess starts much closer to the supersystem solution than SAD and converges in noticeably fewer SCF iterations. For a strongly-coupled / covalently-bonded partition the inter-fragment blocks matter and SAD or SAP is usually better.

Honest errors. The fragments must be a disjoint partition covering every atom exactly once, and their charges must sum to the molecule charge (electron conservation); each violation raises a clear ValueError. A fragment SCF that fails to converge raises; passing fragments= without initial_guess=FRAGMO raises; and FRAGMO on a periodic system raises a clean NotImplementedError. FRAGMO currently requires a named basis set (it rebuilds each fragment’s basis by name). The automatic fragments= assembly belongs to the molecular run_* wrappers. run_neb has no fragments= keyword; molecular NEB accepts FRAGMO only when the caller has already populated complete read_density matrices in the selected SCF options.

(Roadmap D2h, the fragment-MO assembly guess; assembling the fragment densities is equivalent to a block-diagonal MO assembly for seeding and needs no occupation bookkeeping.)

READ: restart from a prior calculation¶

READ skips a from-scratch guess and starts an SCF from a density you already have. It is useful for an explicitly managed geometry sequence, a supported first-band NEB seed, or continuing a job that ran out of iterations. The accepted source depends on the entry point; READ is not a universal warm-start switch for every optimizer and reaction-path driver. Select it with opts.initial_guess = vq.InitialGuess.READ and use one of the sources that the chosen driver documents.

In-memory molecular single point, from a prior result, via the read_from= keyword on the RHF/UHF/RKS/UKS wrappers:

r1 = vq.run_rhf(mol, basis, opts)                     # first point
opts2 = vq.RHFOptions()
opts2.initial_guess = vq.InitialGuess.READ
r2 = vq.run_rhf(mol2, basis, opts2, read_from=r1)     # restart from r1

read_from accepts any prior RHF / UHF / RKS / UKS result. An in-memory result carries no basis or geometry, so this path needs the same basis as the prior run; for a changed basis or geometry, read from a file. run_neb does not accept a read_from= keyword.

From a .qvf or Molden file, via opts.read_path:

opts.read_path = "prev.qvf"      # reads the QVF restart section
# or
opts.read_path = "prev.molden"
# ORCA's default exporter name is accepted directly too:
opts.read_path = "prev.molden.input"

Write the source with QVF output (run_job(..., output_qvf=True)), which embeds wavefunction.gto for molecular jobs and x_vibeqc.bloch_wavefunction for periodic closed-shell jobs, or write a Molden file for molecular / Gamma-style orbital restarts. READ accepts both the raw Gaussian contraction coefficients written by vibe-qc and the primitive-normalized convention written by ORCA’s orca_2mkl. It selects the convention from the MO orthonormality condition and refuses the file if neither interpretation is valid, if an MO omits its occupation or has a non-finite coefficient, or if any occupied MO has an incomplete or duplicate AO-coefficient block, rather than silently changing the electron count. If both contraction interpretations pass, the smaller orthonormality residual wins; an exact tie retains the historical raw-coefficient interpretation. Coefficient tables rounded too coarsely to establish the invariant are rejected rather than guessed. This restart reader currently supports pure-spherical, separate-shell Molden bases; combined sp shells and Cartesian [6D] / [10F] variants are not a supported restart surface.

Closed-shell drivers use the prior total density, open-shell drivers the per-spin pair (opts.read_density / opts.read_density_alpha / opts.read_density_beta).

Projection across a changed basis or geometry. A file source carries its own basis and geometry, so the prior density is projected onto the current basis automatically:

\[ \mathbf{D}_{\rm cur} = \mathbf{P}\,\mathbf{D}_{\rm prior}\,\mathbf{P}^{T}, \qquad \mathbf{P} = \mathbf{S}_{tt}^{-1}\,\mathbf{S}_{tm}, \]

with $\mathbf{S}_{tt}$ the current-basis overlap and $\mathbf{S}_{tm}$ the cross-basis overlap $\langle \chi^{\rm cur} | \chi^{\rm prior} \rangle$. The projection is exact when the source and target bases coincide at the same geometry. It prepares a prior density for the current call’s basis and geometry, so an explicitly managed scan can restart each point from a file. For molecular Gaussian NEB, options.read_path is resolved once against the reactant endpoint to seed the initial band; later outer iterations use NEB’s own warm_start cache. In-memory restarts skip projection and require a matching basis.

Aliases. In the periodic guess-string interface, "moread" and "coread" are accepted ORCA-style spellings of READ; the molecular run_* drivers take the enum directly.

Periodic restart. READ also restarts periodic SCF jobs, including explicit geometry-scan and extrapolation sequences on slabs and bulk. vq.run_periodic_job(..., initial_guess="read", read_from=prior_or_path) takes a prior periodic result (in-memory) or a .qvf / .molden path at Gamma. The prior g=0 cell density is projected onto the current cell basis (the same MINAO-style projector, evaluated with the cell overlaps) and injected at g=0, where it Bloch-sums to a k-independent $\mathbf{D}(\mathbf{k})$ for the first Fock build. Wired on the Ewald, GDF, and BIPOLE Gamma drivers.

Closed-shell multi-k restarts are supported from an in-memory prior result or a current vibe-qc .qvf archive on the GDF, GPW, and GAPW routes. In this case vibe-qc rebuilds every complex Bloch density block from the prior coefficients and occupations, $\mathbf{D}(\mathbf{k}) = \mathbf{C}(\mathbf{k}) f(\mathbf{k}) \mathbf{C}^{\dagger}(\mathbf{k})$, and starts the next multi-k SCF from that per-k density list. Current QVF archives carry the all-k restart data in the first-party vendor section x_vibeqc.bloch_wavefunction; legacy QVF archives with only wavefunction.gto still fail closed because that visualization section carries at most one selected k-point. .molden file sources and open-shell multi-k READ are likewise still pending because they need all-k complex Bloch data, and for open-shell restarts per-spin per-k density blocks.

This run_periodic_job restart surface does not automatically extend to every geometry driver. Fixed-cell BIPOLE optimization rejects READ because it cannot preserve the restart across displaced geometries; run the restart single point first, then optimize from a geometry-defined SAD/HCORE/ATOMSPIN guess. Periodic Gaussian BIPOLE NEB likewise rejects READ and FRAGMO before dry-run: its image and finite-difference displacement loop has neither a geometry-projected full-lattice restart seam nor a fragment-assembly seam, and run_neb has no read_from= or fragments= keyword. Closed-shell RHF/RKS periodic NEB also rejects PATOM because those drivers do not implement the in-field seed; UHF/UKS may use PATOM or ATOMSPIN. Run any requested restart or fragment guess as a supported single point first, then begin the NEB band from SAD/HCORE (or an unrestricted magnetic seed).

Like HCORE, READ carries no citation: it is mechanical, not a published method.

When to use which, decision table¶

Situation	Recommended	Why
Default, any molecular system	AUTO	Engine picks SAP for closed-shell light, SAD otherwise
Default, any periodic system	AUTO	Resolves to SAD; fixes ionic-insulator failure modes
Difficult open-shell convergence	SAD	Robust; OH/6-31G* false-minimum example
Heavy atoms (Z > 36), closed shell	SAP	Cheaper than SAD’s atomic-SCF cost
Transition-metal complex	SAD (AUTO default) or PATOM	PATOM re-polarises the SAD densities, helping the d-shell ordering that SAP’s closed-shell averaging mis-orders
Unit tests, deterministic reference	HCORE	System-independent, no implementation-dependent atomic SCF
Explicit molecular/periodic scan sequence	READ	Restart each supported single-point call from the previous result or a projected file source
Molecular Gaussian NEB first-band seed	READ via `options.read_path` or complete preloaded densities	`run_neb` has no `read_from=` keyword; its internal `warm_start` handles later outer iterations
BIPOLE fixed-cell optimization	SAD/HCORE/ATOMSPIN, not READ	The optimizer rejects READ across displaced geometries; run a READ single point first if needed
Weakly-interacting dimer / vdW complex / cluster	FRAGMO	Superpose converged fragment densities; starts close to the supersystem solution, fewer iterations than SAD
Broken-symmetry singlet	SAD + `guessmix_angle_deg` (when it ships)	α HOMO/LUMO rotation to break spatial symmetry

The remaining method marked “when it ships” (broken-symmetry guessmix_angle_deg) is still scaffolded in the engine; its concrete implementation is tracked in docs/roadmap.md under the initial-guess track (§G2) and the DIIS-family section (§D1).

ROHF and ROKS: a smaller supported set¶

The ROHF / ROKS drivers are the pure-Python Roothaan loop, which seeds from per-spin densities only. The overlap-dependent guesses (SAP, PATOM, HUECKEL, MINAO) need S inside the guess engine and are therefore not available there. ROHFOptions.initial_guess and ROKSOptions.initial_guess accept:

Value	Effect
`"sad"` / `InitialGuess.SAD`	Superposition of atomic densities (default)
`"hcore"` / `"core"` / `InitialGuess.HCORE`	Core-Hamiltonian guess
`"auto"` / `InitialGuess.AUTO`	Resolves to SAD for an open-shell system

Both spellings (the InitialGuess enum and the canonical string) are accepted and normalised internally. Anything else raises ValueError:

from vibeqc import InitialGuess
from vibeqc.rohf import ROHFOptions, run_rohf

opts = ROHFOptions()
opts.initial_guess = InitialGuess.HCORE   # enum: honoured
opts.initial_guess = "core"               # string: same guess

opts.initial_guess = InitialGuess.SAP     # ValueError: not supported
opts.initial_guess = "bogus"              # ValueError: unknown guess

This is a behaviour change. Previously any value other than the string "core" was silently ignored here and the driver used SAD regardless, so assigning the InitialGuess enum looked like it took effect but did not. A pinned guess is now either honoured or rejected, never silently substituted.

The SAD spin density is rescaled to the target $S_z$ here. The Hund-split atomic densities described under Multiplicity and the initial guess carry the free-atom spin populations, which for a molecule need not sum to $n_\alpha - n_\beta$: NH$_2(^2B_1)$ and CH$_3$ get 7α/2β, OH 6α/3β, O$_2(^3\Sigma)$ 10α/6β. The unrestricted drivers refill by aufbau on their first cycle, so it never reaches their fixed point. The Roothaan loop cannot: it builds its closed / open / virtual projectors from the density itself, and hands the same iterate to the accelerator history. A density-mode ROHF / ROKS guess therefore keeps its total density and Hund spin pattern but has that pattern’s magnitude scaled so $\mathrm{Tr}(S(D_\alpha - D_\beta)) = n_\alpha - n_\beta$. Before that, NH$_2$/STO-3G ROHF on plain defaults converged, reporting success with no warning, to the excited $^2A_1$ state 95 mHa above the ground state (issue #119).

Worked examples¶

Comparing iteration counts across guesses¶

import vibeqc as vq
from vibeqc import Atom, Molecule, BasisSet, RHFOptions, run_rhf

mol = Molecule([
    Atom(8, [0.0,  0.0,  0.0]),
    Atom(1, [0.0,  1.43, -0.98]),
    Atom(1, [0.0, -1.43, -0.98]),
])

for basis_name in ["sto-3g", "6-31g*", "cc-pvdz"]:
    basis = BasisSet(mol, basis_name)
    for kind in [vq.InitialGuess.HCORE,
                 vq.InitialGuess.SAD,
                 vq.InitialGuess.SAP]:
        opts = RHFOptions()
        opts.initial_guess = kind
        opts.conv_tol_energy = 1e-10
        r = run_rhf(mol, basis, opts)
        print(f"{basis_name:>10}  {kind.name:<6}  "
              f"n_iter={r.n_iter:2d}  E={r.energy:.6f}")

Output on H₂O (typical):

    sto-3g  HCORE   n_iter= 8  E=-74.964156
    sto-3g  SAD     n_iter= 8  E=-74.964156
    sto-3g  SAP     n_iter= 8  E=-74.964156
    6-31g*  HCORE   n_iter=12  E=-76.006529
    6-31g*  SAD     n_iter=10  E=-76.006529
    6-31g*  SAP     n_iter=11  E=-76.006529
   cc-pvdz  HCORE   n_iter=12  E=-76.023841
   cc-pvdz  SAD     n_iter=11  E=-76.023841
   cc-pvdz  SAP     n_iter=11  E=-76.023841

All three converge to the same energy to machine precision; iteration counts shift by ±1-2 depending on basis. SAP’s advantage grows with atom heaviness, Lehtola’s Tables I-II show 40-60% iteration-count reductions on Au / Ag / I-containing systems.

Open-shell radical (OH)¶

The textbook example where HCore gets the wrong minimum:

from vibeqc import Atom, Molecule, BasisSet, UHFOptions, run_uhf
import vibeqc as vq

A2B = 1.0 / 0.529177210903    # Angstrom → bohr
mol = Molecule(
    [Atom(8, [0.0, 0.0, 0.0]),
     Atom(1, [0.97 * A2B, 0.0, 0.0])],
    multiplicity=2,
)
basis = BasisSet(mol, "6-31g*")

# HCore lands on a false minimum 0.16 Ha above the true one:
opts = UHFOptions()
opts.initial_guess = vq.InitialGuess.HCORE
r_hcore = run_uhf(mol, basis, opts)
print(f"HCORE: E = {r_hcore.energy:.6f}  S² = {r_hcore.s_squared:.3f}")
# HCORE: E = -75.221182  S² = ~0.764 (or worse)

# SAD finds the right minimum, S² near 0.755 (doublet):
opts.initial_guess = vq.InitialGuess.SAD
r_sad = run_uhf(mol, basis, opts)
print(f"SAD:   E = {r_sad.energy:.6f}  S² = {r_sad.s_squared:.3f}")
# SAD:   E = -75.380931  S² = 0.755

run_uhf defaults to AUTO, which resolves to SAD for any open-shell input, so omitting opts.initial_guess entirely also reaches the correct minimum.

Periodic ionic insulator (NaCl)¶

import vibeqc as vq
import numpy as np
from vibeqc import Atom, BasisSet, PeriodicSystem, PeriodicRHFOptions

a = 10.7      # NaCl rocksalt lattice parameter, bohr (close enough for sto-3g)
lat = a * np.eye(3)
atoms = [
    Atom(11, [0.0, 0.0, 0.0]),         # Na
    Atom(17, [a/2, a/2, a/2]),         # Cl
]
sysp = PeriodicSystem(3, lat, atoms, charge=0, multiplicity=1)
basis = BasisSet(sysp.unit_cell_molecule(), "sto-3g")

# AUTO → SAD for periodic. The default works.
opts = PeriodicRHFOptions()
opts.conv_tol_energy = 1e-8
result = vq.run_rhf_periodic_gamma(sysp, basis, opts)
print(f"NaCl Γ-only RHF: E = {result.energy:.6f}, "
      f"n_iter = {result.n_iter}")

# If you explicitly request HCORE, you may hit the bombing failure
# mode (first-iter energy ~+30 000 Ha, recovers slowly or not at all):
opts.initial_guess = vq.InitialGuess.HCORE
# ... (don't do this for ionic insulators)

The SAD default closes the v0.5.6 NaCl-bombing failure mode that prompted this whole refactor.

Inspecting what AUTO chose¶

The resolved kind appears in the SCF log:

[scf]  initial guess: SAP (Lehtola/Visscher/Engel 2020, sap_helfem_large)

To check programmatically without running SCF, call the engine directly through its private pybind binding:

from vibeqc._vibeqc_core import _guess_closed_shell_density, InitialGuess

# Pass kind=AUTO; the engine returns the density for closed-shell
# or None for HCORE - either way you can tell from the SCF log
# what kind was resolved.
D = _guess_closed_shell_density(
    mol, basis, mol.n_electrons() // 2,
    InitialGuess.AUTO, is_periodic=False,
)

This is mainly useful for debugging / scripting; production code should just set opts.initial_guess = vq.InitialGuess.AUTO and read the log.

Diagnosing convergence problems¶

When SCF won’t converge or lands on the wrong state, the initial guess is one of three places to look (the others are damping/DIIS and the system geometry). Quick checklist:

Is the energy diverging or oscillating wildly? Switch to SAD. If you’re already on SAD or AUTO and it’s still bad, the problem isn’t the guess, check the geometry and damping (scf_convergence.md).
Is it converging but to the wrong S²? Open-shell trap. Try SAD (UHF default) and check result.s_squared against the expected $S(S+1)$ from the multiplicity. For broken-symmetry singlets, use SAD and (once it ships) guess_request.guessmix_angle_deg = 45.
Iteration count creeping up over a geometry scan? Use READ to restart from the previous step’s density (read_from= the prior result, or read_path= a .qvf); it projects across the geometry change. See the READ section above.
Periodic system with a charge density that looks “ionic” but the metal sites are wrong? SAD is doing its job, but the AUTO- selected even-spin split may not bias correctly for AFM / FM states. Use SAD explicitly + (future) atomic_occupations to set per-element starting occupations.

Extending the engine¶

Adding a new method is a small follow-up commit on the existing scaffold. The contract:

Declare it. Add the enum value in cpp/include/vibeqc/guess.hpp (the nine current kinds are all implemented; this step is for a future tenth).
Implement it. Add a free function in cpp/src/guess.cpp (or a new file cpp/src/guess_<name>.cpp and add it to the CMakeLists) that builds whichever artifact the method produces (density, Fock, or MOs).
Wire it. Add a case InitialGuess::<NAME>: arm in GuessEngine::build_closed_shell (and build_open_shell if it handles spin); populate the matching field on GuessResult and write a provenance string.
Test it. Add at least one parity test (tests/test_guess.py): the new method must converge to the same total energy as HCORE / SAD on a small system to machine precision.

See the existing SAP implementation (commit d2fbfe9) for a worked example. The whole landing was ~500 lines including tests, citation entry in docs/license.md, and CHANGELOG.

Result types, closed-shell vs open-shell¶

The engine returns one of two tagged structs depending on which entry point you call. Each result has at most one of the three “modes” populated; the rest are empty.

`GuessClosedShellResult`¶

struct GuessClosedShellResult {
    InitialGuess resolved_kind;  // what AUTO actually picked
    std::string provenance;      // human-readable trace
    Eigen::MatrixXd D;           // density-mode (SAD, PATOM, HUECKEL, MINAO)
    Eigen::MatrixXd F_guess;     // Fock-mode (HCORE, SAP)
    Eigen::MatrixXd C_guess;     // MO-mode (reserved; no current guess)
    Eigen::VectorXd eps_guess;
};

`GuessOpenShellResult`¶

struct GuessOpenShellResult {
    InitialGuess resolved_kind;
    std::string provenance;
    Eigen::MatrixXd D_alpha, D_beta;               // density-mode
    Eigen::MatrixXd F_guess_alpha, F_guess_beta;   // Fock-mode
    Eigen::MatrixXd C_guess_alpha, C_guess_beta;   // MO-mode
    Eigen::VectorXd eps_guess_alpha, eps_guess_beta;
};

The Python helpers in vibeqc.guess (initial_density_closed_shell, initial_densities_open_shell) already unpack the result for the Python periodic SCF drivers; most users won’t touch the result type directly.

`InitialGuess` enum¶

enum class InitialGuess {
    AUTO,      // default: engine picks per system
    HCORE,     // F = T + V_ne
    SAD,       // superposition of atomic densities
    SAP,       // superposition of atomic potentials (Lehtola 2020)
    PATOM,     // SAD + one in-field re-polarisation step
    HUECKEL,   // parameter-free generalised Wolfsberg-Helmholz
    MINAO,     // minimal-AO (ANO-RCC) reference projection
    READ,      // restart from a prior result / .qvf / .molden
    FRAGMO,    // superposition of converged fragment densities (molecular)
};

Bound to Python as vibeqc.InitialGuess.

References¶

S. Lehtola, L. Visscher, E. Engel, Efficient Implementation of the Superposition of Atomic Potentials Initial Guess for Electronic Structure Calculations in Gaussian Basis Sets, J. Chem. Phys. 152, 144105 (2020). DOI:10.1063/5.0004046
S. Lehtola, Assessment of Initial Guesses for Self-Consistent Field Calculations, J. Chem. Theory Comput. 15, 1593 (2019). DOI:10.1021/acs.jctc.8b01089
J. H. Van Lenthe, R. Zwaans, H. J. J. Van Dam, M. F. Guest, Starting SCF Calculations by Superposition of Atomic Densities, J. Comput. Chem. 27, 926 (2006). DOI:10.1002/jcc.20393
Q. Sun et al., PySCF: the Python-based simulations of chemistry framework, WIREs Comput. Mol. Sci. 8, e1340 (2018), description of the minao projection used in PySCF.