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      # (planned) SAD + in-field re-polarisation
vq.InitialGuess.HUECKEL    # (planned) extended Hückel
vq.InitialGuess.MINAO      # (planned) minimal-AO projection
vq.InitialGuess.READ       # (planned) read from a prior calculation

Three are implemented in v0.9 (HCORE, SAD, SAP). The remaining four raise a clear std::runtime_error pointing at the roadmap entry — they are scaffolded in the engine so adding them is a small follow-up per method.

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`	`D`
Fock	`HCORE`, `SAP`, `HUECKEL`	`F_guess`
MO	`MINAO`, `READ`	`C_guess`, `eps_guess`

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.

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 Aufbau occupations of the atom’s orbitals (degenerate orbitals receive equal fractional occupations so the result is spherical). 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; PATOM will be the right answer once it ships
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 …)"

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 (today); PATOM (when it ships)	SAP’s closed-shell averaging mis-orders d-shell occupations
Unit tests, deterministic reference	HCORE	System-independent, no implementation-dependent atomic SCF
Geometry scan / NEB / opt	READ (when it ships)	Restart from the previous step’s orbitals
Broken-symmetry singlet	SAD + `guessmix_angle_deg` (when it ships)	α HOMO ↔ LUMO rotation to break spatial symmetry

The “today” column is what the v0.9 release ships. Methods marked “when it ships” are scaffolded in the engine; their concrete implementations are tracked in docs/roadmap.md under “v0.6.x — initial-guess overhaul” (§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:

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 (once it ships) to restart from the previous step’s MOs. Until then, expect 8–15 iterations per point.
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 (already present for the planned methods — PATOM, HUECKEL, MINAO, READ).
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)
    Eigen::MatrixXd F_guess;     // Fock-mode (HCORE, SAP, HUECKEL)
    Eigen::MatrixXd C_guess;     // MO-mode (MINAO, READ)
    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,     // (planned) SAD + in-field re-polarisation
    HUECKEL,   // (planned) extended Hückel minimal basis
    MINAO,     // (planned) minimal-AO projection
    READ,      // (planned) restart from a prior calc
};

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.