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 Γ-only 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 eight 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

All eight are implemented today. READ works for molecular and Γ-point periodic SCF (a multi-k READ raises — the QVF wavefunction.gto section is Γ-only; see the READ section). PATOM, HUECKEL, and MINAO are molecular-only: a periodic SCF selecting one raises a clear error pointing at the roadmap (“not yet implemented for periodic”).

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, HUECKEL, MINAO

D

Fock

HCORE, SAP

F_guess

READ sits outside this table: it is resolved in the run_* wrapper rather than the engine, injecting a prior density as opts.read_density (closed shell) or opts.read_density_alpha / opts.read_density_beta (open shell). The READ section below covers it.

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, spatially symmetric. Given \(n_\alpha\) and \(n_\beta\), the open-shell guess produces a starting point with the right counts but the same spatial shape for both spins:

  • SAD (the AUTO choice for any open shell) uses a proportional split, \(\mathbf{D}_\alpha = (n_\alpha/n_e)\,\mathbf{D}_{\rm SAD}\) and \(\mathbf{D}_\beta = (n_\beta/n_e)\,\mathbf{D}_{\rm SAD}\).

  • The UKS packaging diagonalises \(\mathbf{F}_{\rm SAD}\) once and keeps the first \(n_\alpha\) and \(n_\beta\) columns: the same orbitals filled to two different depths.

Either way the starting alpha and beta densities are spatially identical. The physical spin density develops during SCF: once \(n_\alpha \ne n_\beta\) the per-spin Fock matrices differ (\(\mathbf{F}_\alpha \ne \mathbf{F}_\beta\)) and the asymmetry grows iteration by iteration. 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 symmetric start biases toward the most symmetric one.

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 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 drivers (a two-phase SCF that locks a high-spin state, then restarts at the target multiplicity from that density via the READ machinery). Periodic SPINLOCK on a driver that does not implement the requested mode (GDF, BIPOLE, the multi-k UHF Ewald path, or SPIN_SCHEDULE on multi-k UKS) fails closed with a clear error rather than silently ignoring the request.

Still pending: BETALOCK (locking only the beta count), SPIN_SCHEDULE on the multi-k and GDF/BIPOLE 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 not yet implemented. A periodic SCF with opts.initial_guess = SAP raises a clear error pointing at the roadmap. Lattice-summed erf_nuclear integrals are needed; the work slots in alongside the existing Γ-only periodic JKBuilder interface.

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

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

Periodic SAP not yet wired; SAD fixes the NaCl-bombing failure mode

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 (molecular)

Three further density-mode guesses landed for molecules; each returns a starting density the SCF consumes directly, and each raises a clear “not yet implemented for periodic” error on a periodic system.

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

READ: restart from a prior calculation

READ skips a from-scratch guess and starts the SCF from a density you already have. It is the right choice for a geometry optimisation, an NEB chain, a relaxed scan, or continuing a job that ran out of iterations: each step starts from the last converged density and needs only a few iterations. Select it with opts.initial_guess = vq.InitialGuess.READ and give it one of three sources.

In-memory, from a prior result, via the read_from= keyword:

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.

From a .qvf or .molden, via opts.read_path:

opts.read_path = "prev.qvf"      # reads the wavefunction.gto section
# or
opts.read_path = "prev.molden"

Write the source with QVF output (run_job(..., output_qvf=True), which embeds wavefunction.gto) or a Molden file so the basis and MO coefficients are there to restart from.

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 bases coincide, and is what lets READ carry a density across the geometry steps of an optimisation or NEB, or across a basis change. 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 (Γ-point) restart. READ also restarts a Γ-point periodic SCF — the high-value NEB / geometry-scan / extrapolation case 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; 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 Γ drivers. Multi-k per-k restart needs complex Bloch coefficients (the QVF wavefunction.gto section is Γ-only) and is out of scope — a multi-k READ raises with a roadmap pointer.

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

Geometry scan / NEB / opt

READ

Restart from the previous step’s orbitals (projected across the geometry change)

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

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:

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

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

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

  4. 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:

  1. Declare it. Add the enum value in cpp/include/vibeqc/guess.hpp (the eight current kinds are all implemented; this step is for a future ninth).

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

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

  4. 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
};

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.

See also