Initial Guess

Every SCF (HF/DFT) calculation needs a starting set of molecular orbitals — the initial guess. A good guess can be the difference between converging in a handful of iterations and not converging at all, and for excited-state or broken-symmetry work the guess often determines which solution you find. The guess strategy is selected with the guess keyword.

guess <strategy> [extra arguments]

Available strategies

`guess …`	Strategy	Notes
(omitted) or `generate`	Maximum-overlap MO guess (default)	Per-element guess orbitals from the basis-set library; no extra files needed.
`<file>` / `<ca> <cb>`	Read orbitals from file	Binary `c0`/`ca0`/`cb0` orbital files, or a `.molden` file.
`hcore`	Core-Hamiltonian guess	Diagonalize the one-electron (core) Hamiltonian. Cheap, usually a poor starting point.
`sad [file]`	SAD — densities computed on the fly	Superposition of Atomic Densities; runs a small atomic SCF per unique element. Works with any basis.
`sadlp [file]`	SAD — precomputed densities from disk	Loads stored atomic densities shipped with TeraChem. Limited basis sets/elements (see caveat).
`frag <file>`	Fragment guess	Superposition of converged SCF densities of user- or auto-defined molecular fragments.
`project`	Wavefunction projection	Project a converged small-basis solution into the target basis. See Wavefunction Projection.
`sap`	SAP — superposition of atomic potentials	Builds an effective potential from tabulated atomic potentials and diagonalizes it.
`exciton <file>`	Exciton model	Fragment-localized guess for the ab initio exciton model. See Exciton Model.

The default is not SAD

With no guess keyword (or guess generate) TeraChem uses a maximum-overlap molecular-orbital guess built from per-element guess orbitals stored alongside each basis set. If you are used to codes that default to SAD, request guess sad explicitly.

The default guess: maximum-overlap atomic orbitals (`generate`)

With no guess keyword (or guess generate), TeraChem builds the starting orbitals directly from precomputed atomic orbitals — without any molecular SCF iterations or Fock-matrix diagonalization. The algorithm is the Hartree–Fock part (denoted HF-INIT) of the GVB-INIT scheme of Langlois, Yamasaki, Muller, and Goddard.¹

Only the HF part of GVB-INIT is used

The reference describes a full initial-guess procedure for generalized valence bond (GVB) wave functions, of which HF-INIT is the first stage. TeraChem uses only HF-INIT — it constructs the occupied HF molecular orbitals and stops. The subsequent GVB steps (piecewise bond/lone-pair localization and the second natural orbitals) are not performed.

Stored atomic orbitals

For each element, the basis-set library (the basis.ao files; see Basis Sets) stores a fixed set of atomic guess orbitals, partitioned into core orbitals $\{\theta^{\text{core}}_i\}$ and valence orbitals $\{\theta^{\text{val}}_i\}$, each a fixed contraction of that atom's contracted basis functions $\{\chi_\mu\}$. Placing these atom-centred orbitals on every atom and expressing them in the molecular AO basis gives a rectangular coefficient matrix $\mathbf{T}$:

\[ \theta_i(\mathbf{r}) = \sum_\mu T_{\mu i}\,\chi_\mu(\mathbf{r}). \]

Valence orbitals — select the maximally-overlapping combinations

The central idea is that the occupied valence molecular orbitals can be chosen as the combinations of valence atomic orbitals with the largest mutual overlap, since these are the incipient bonding orbitals. TeraChem builds the overlap matrix of the valence atomic orbitals in the molecule,

\[ S^{\text{val}}_{ij} = \langle \theta^{\text{val}}_i \,|\, \theta^{\text{val}}_j\rangle = \bigl(\mathbf{T}^{\text{val}\,\mathsf{T}}\,\mathbf{S}\,\mathbf{T}^{\text{val}}\bigr)_{ij}, \]

where $\mathbf{S}$, with $S_{\mu\nu}=\langle\chi_\mu|\chi_\nu\rangle$, is the AO overlap matrix of the whole molecule. This is Eq. (3.1) of Ref. ¹. Diagonalizing,

\[ \mathbf{S}^{\text{val}} = \mathbf{U}\,\mathbf{s}\,\mathbf{U}^{\mathsf{T}}, \qquad \phi^{\text{val}}_j = \frac{1}{\sqrt{s_j}}\sum_i \theta^{\text{val}}_i\,U_{ij}, \]

gives orthonormal valence MOs $\phi^{\text{val}}_j$. The eigenvalue $s_j$ measures how strongly the atomic orbitals reinforce in combination $j$: a large $s_j$ marks a strongly-overlapping (bonding) combination. TeraChem keeps the $N_{\text{occ}}-N_{\text{core}}$ valence MOs with the largest eigenvalues, where $N_{\text{occ}}$ is the number of occupied orbitals (fixed by the electron count and charge) and $N_{\text{core}}$ the number of stored core orbitals.

Overlap, not Fock

Diagonalizing the overlap matrix — rather than a Fock or kinetic-energy matrix as in extended-Hückel-type guesses — is the distinctive feature of this method; in TeraChem's testing it outperforms those alternatives.

Core orbitals

Each stored core atomic orbital is taken directly as a doubly-occupied guess MO (step γ1 of Ref. ¹); no selection is applied.

Assembly and orthonormalization

The $N_{\text{core}}$ core orbitals and the selected valence orbitals are stacked (core first) into a combined set with AO-coefficient matrix $\mathbf{A}$. Because the two sets are not mutually orthogonal, the combined set is orthonormalized through the eigendecomposition of its own overlap matrix, $\tilde{\mathbf{S}}=\mathbf{A}^{\mathsf{T}}\mathbf{S}\,\mathbf{A} = \mathbf{V}\,\boldsymbol{\sigma}\,\mathbf{V}^{\mathsf{T}}$:

\[ \mathbf{C} = \mathbf{A}\,\mathbf{V}\,\boldsymbol{\sigma}^{-1/2}, \qquad \mathbf{C}^{\mathsf{T}}\mathbf{S}\,\mathbf{C} = \mathbf{I}, \]

yielding the final occupied guess-MO coefficients $\mathbf{C}$. The restricted guess density follows as $P_{\mu\nu} = 2\sum_{p}^{\text{occ}} C_{\mu p}C_{\nu p}$; for unrestricted methods the whole procedure is run separately with $N_\alpha$ and $N_\beta$ occupied orbitals.

Relationship to SAD

This guess is closely related to SAD — both assemble the molecule from atomic information without a molecular SCF — but it is not identical. SAD builds the guess from a superposition of atomic densities; the generate guess instead selects molecular orbitals by maximal overlap of the stored atomic orbitals (HF-INIT), which Ref. ¹ explicitly distinguishes from the density-superposition variant (their "HF-DENS"). For typical organic molecules both reproduce the converged SCF occupied orbitals to better than ~98%.

Cartesian basis only

Like sadlp, the generate guess requires a Cartesian basis; it aborts on a true-spherical-basis (true_spherical_basis) calculation.

Reading orbitals from a file

At the end of every SCF run TeraChem writes the converged orbitals to the scratch directory (scr/c0, or scr/ca0 and scr/cb0 for unrestricted methods). Point guess at such a file to restart from it:

guess           scr/c0                # restricted
guess           scr/ca0   scr/cb0     # unrestricted: alpha then beta

If only the alpha file is given for an unrestricted calculation, the beta orbitals are copied from alpha. A .molden file may be supplied instead of a binary orbital file; TeraChem detects the .molden extension and reorders the orbitals automatically.

SAD — superposition of atomic densities

SAD builds the guess density as a direct sum of spherically-averaged atomic densities. It is robust and basis-agnostic, and is the recommended general choice. TeraChem provides two SAD implementations.

`guess sad` — atomic densities computed on the fly

For each unique atom (same element, charge, and spin) TeraChem runs a small atomic UHF/UKS SCF and reuses the result for all equivalent atoms. Because the densities are generated at run time, this variant works with any basis set.

Keyword	Type	Default	Description
`sad_fxnl`	string	`hf`	Functional used for the atomic densities. `hf` or `blyp`.
`sad_spherical`	bool	`yes`	Spherically symmetrize the atomic densities.
`sad_avgspin`	bool	`yes`	Average the alpha and beta atomic densities.

By default each atom is taken neutral, spin-up, with the NIST ground-state spin multiplicity (tabulated up to Z = 86). To override charges and spins per atom, supply a file as the third token. The file has one line per atom, in coordinate-file order, each line giving net charge, spin multiplicity, and spin direction (a = alpha/up, b = beta/down):

spinfile (for CO₂, atoms O C O)

0  3  a
0  3  a
0  1  a

sad.in

coordinates   CO2.xyz
basis         3-21g
method        uhf
charge        0
spinmult      1

guess         sad   spinfile
sad_spherical yes
sad_avgspin   yes

run           energy
end

`guess sadlp` — precomputed atomic densities from disk

sadlp skips the atomic SCF and instead loads tabulated spherical atomic densities shipped with TeraChem from $TeraChem/rhosph/<basis>/<element>.p. This is faster but only available for the basis sets and elements that were pre-tabulated.

Limited and unverified coverage

sadlp ships precomputed UHF densities for only a handful of basis sets:

3-21G and 6-31G — H through Zn
6-31G*, 6-31G**, 6-31+G*, 6-31++G** — H through Ca

It additionally requires a Cartesian basis (it aborts on true_spherical_basis) and that the $TeraChem environment variable points at an install containing the rhosph/ density library. Unlike the other strategies, sadlp is not exercised by the integration-test suite, so treat it as a legacy convenience and prefer guess sad for new work.

SAP — superposition of atomic potentials

guess sap constructs an effective one-electron potential from tabulated atomic potentials and diagonalizes it. The atomic potentials come from a dedicated basis-like library, selected with sap_basis (default sap_grasp_large, based on GRASP relativistic atomic data).

Keyword	Type	Default	Description
`sap_basis`	string	`sap_grasp_large`	Name of (or path to) the SAP potential library.

sap.in

coordinates   CO2.xyz
basis         6-31gs
method        uhf
charge        0
spinmult      1

guess         sap
sap_basis     sap_grasp_large

run           energy
end

Fragment guess

guess frag partitions the system into molecular fragments, converges an independent (UHF-like) SCF on each fragment, and superimposes the resulting fragment densities. This is especially useful for getting broken-symmetry or charge/spin-localized starting points that SAD cannot represent — e.g. forcing a particular fragment to carry a given charge or spin.

The fragments can be specified explicitly in a file, or detected automatically.

Defining fragments from a file

Supply the fragment file as the third token (guess frag <file>). The first line is the number of fragments; each subsequent line describes one fragment with four tokens — number of atoms, net charge, spin multiplicity, and spin direction (a/b) — with atoms assigned to fragments in coordinate-file order:

fragfile (water dimer, two 3-atom fragments)

2
3 0 1 a
3 0 1 a

The a/b spin-direction label sets the sign of each fragment's $S_z$, so a pair such as … a / … b produces an overall low-spin (antiferromagnetically coupled) starting guess. For ethane split into two methyl radicals coupled to a singlet:

ethane two CH₃ fragments

2
4 0 2 a
4 0 2 b

frag.in

coordinates   water2.xyz
basis         3-21g
method        rhf
charge        0
spinmult      1

guess         frag   fragfile

run           energy
end

Automatic fragment detection (`frag_auto`)

Instead of a file, TeraChem can determine the fragments itself. frag_auto selects the method:

`frag_auto`	Fragmentation method
`amino`	Detect amino-acid residues (and waters) from a PDB coordinate file. Charges are set from protonation state; spins are minimized; non-terminal residues are assigned triplet multiplicity to push unpaired electrons to opposite ends.
`kmeans` (or `true`)	Spatial k-means clustering. Charges assumed zero; spins minimized.
`spec`	Spectral clustering. Charges assumed zero; spins minimized.

frag_auto with amino-acid detection

coordinates   protein.pdb
basis         3-21g
method        rhf
charge        0

guess         frag
frag_auto     amino

run           energy
end

frag_auto with k-means clustering

coordinates       water2.xyz
basis             3-21g
method            rhf

guess             frag
frag_auto         kmeans
frag_auto_numfrags 2

run               energy
end

Fragment-guess keywords

Keyword	Type	Default	Description
`frag_auto`	string	(off)	Auto-detect fragments: `amino`, `kmeans`, or `spec`. Disregards the fragment file.
`frag_auto_numfrags`	int	$\sqrt{N_\text{atoms}}$	Target number of fragments for `kmeans`/`spec` (fewer may result if empty clusters appear).
`frag_guess`	string	`generate`	Initial guess used within each fragment SCF: `generate` or `hcore`.
`frag_minspin`	bool	`no`	Ignore the multiplicities in the fragment file and minimize each fragment's multiplicity from NIST ground-state spins.
`frag_fxnl`	string	`hf`	Functional for the fragment densities: `hf` or `blyp` (mirrors `sad_fxnl`).
`frag_avgspin`	bool	`yes`	Average alpha/beta fragment densities (mirrors `sad_avgspin`).
`frag_dblprec`	bool	`no`	Use double precision when converging the fragment SCFs.
`frag_maxit`	int	`50`	Maximum SCF iterations per fragment.
`frag_threshscale`	float	`1.0`	Scale factor applied to the SCF convergence threshold for fragments.

Embedding fragments in an MM-charge field (`frag_mm`)

The fragment densities can be refined self-consistently in the electrostatic field of the surrounding fragments: after an initial vacuum calculation, point charges derived from each fragment are used as a field in which the fragments are recomputed.

Keyword	Type	Default	Description
`frag_mm`	bool	`no`	Recompute fragment densities in the field of charges derived from the other fragments.
`frag_mm_chargetype`	string	`vdd`	Charge model for the embedding field: `vdd` (`1`) or `mulliken` (`0`).
`frag_mm_coeff`	float	`1.0`	Scale factor applied to the embedding charges.
`frag_mm_multipass`	int / `converge`	`2`	Number of charge/density refinement passes; `converge` iterates until every fragment needs ≤ 3 SCF iterations.

Breaking spin symmetry for unrestricted (broken-symmetry) solutions

For an open-shell singlet — or any system whose lowest UHF/UKS solution is spin-broken (e.g. an antiferromagnetically-coupled diradical, a stretched bond, or a transition-metal complex) — the symmetric guess is often a saddle point: with identical α and β orbitals the SCF stays trapped at the restricted (RHF/RKS) solution and never finds the lower broken-symmetry one. TeraChem offers two ways to perturb the guess so the SCF can drop into the broken-symmetry minimum.

`mixguess` — rotate the α/β HOMO and LUMO

mixguess mixes the HOMO and LUMO of the unrestricted guess, rotating the α pair by +mixguess degrees and the β pair by −mixguess degrees. Because the two spin channels are rotated in opposite directions, the resulting α and β orbitals are no longer spatially identical, which lets the SCF break symmetry.

Keyword	Type	Default	Description
`mixguess`	float (degrees)	`45` for singlets, `0` otherwise	Rotation angle between the HOMO and LUMO in the unrestricted initial guess. The α orbitals are rotated by `+mixguess`, the β orbitals by `−mixguess`. Set `0` to disable.

Only meaningful for unrestricted methods

HOMO/LUMO mixing only makes sense for unrestricted (UHF/UKS) calculations, where α and β orbitals are independent. For singlets TeraChem enables it by default (45°), precisely because the symmetric singlet guess is the case most prone to collapsing onto the restricted solution; for non-singlets it is off by default. It has no effect on restricted (RHF/RKS) runs.

MOM turns mixing off

Maximum-overlap-method calculations (mom/imom) disable HOMO/LUMO guess mixing automatically, so that the occupation requested in the $excitations block is not scrambled. See MOM / IMOM.

Asymmetric level shifting

Level shifting (levelshift yes) can also be used to favor a broken-symmetry solution: by applying different shifts to the α and β virtual manifolds with levelshiftvala and levelshiftvalb, the two spin channels see different effective Fock operators and are pushed toward distinct spatial orbitals — a truly unrestricted (spin-polarized) solution — rather than relaxing back to the symmetric one.

UHF singlet, mixing + asymmetric level shift

coordinates     diradical.xyz
basis           6-31gss
method          uhf
charge          0
spinmult        1

guess           sad
mixguess        45            # rotate alpha/beta HOMO-LUMO by ±45°
levelshift      yes
levelshiftvala  0.3           # different alpha ...
levelshiftvalb  0.1           # ... and beta shifts break the symmetry

run             energy
end

See the SCF Overview for the full description of the level-shift keywords.

Practical notes and gotchas

Semiempirical / xTB override the guess

When running GFN-xTB, GFN2-xTB, or other semiempirical methods, any guess other than read-from-file is silently replaced by the Hcore/extended-Hückel guess. A guess sad requested alongside xtb will not take effect.

Projection disables guess purification

guess project automatically turns off density-matrix purification of the guess (the two are currently incompatible). This is expected behavior, not an error.

The default purify guess step (purifying the initial density before the SCF) applies to the density-based guesses (sad, sadlp, frag); set purify no to disable it, as in the worked examples above.

J.-M. Langlois, T. Yamasaki, R. P. Muller, W. A. Goddard III, "Rule-Based Trial Wave Functions for Generalized Valence Bond Theory," J. Phys. Chem. 98, 13498–13505 (1994). TeraChem implements the Hartree–Fock initial-guess (HF-INIT) portion, §3.1. ↩↩↩↩