Skip to content

SCF convergence, initial guesses, and TRAH: getting the reference to converge

Every wavefunction and DFT calculation in OpenQP begins by solving the self-consistent field (SCF) equations: you guess a set of molecular orbitals, build the Fock matrix they imply, diagonalize it to get new orbitals, and repeat until the orbitals stop changing. Most closed-shell molecules converge in a dozen iterations without any thought. But stretched bonds, transition-metal complexes, diradicals, and near-degenerate states can make the default extrapolation stall, oscillate, or slide into the wrong solution. This tutorial shows the two knobs that fix almost all of those cases — the initial guess ([guess] type) and the converger ([scf] converger_type, including SOSCF and TRAH) — on a plain water/PBE calculation you can run in a second.

A little theory

The SCF cycle is a fixed-point iteration, and a naive "plug the new orbitals straight back in" scheme frequently diverges. OpenQP gives you a ladder of convergers, each doing more work per iteration but converging more robustly:

  • DIIS (the default) extrapolates the next Fock matrix from a history of previous ones. Cheap and excellent for well-behaved, closed-shell cases.
  • SOSCF (second-order SCF) uses curvature — an approximate orbital-rotation Hessian — to take Newton-like steps once you are near the solution, so it finishes the tail of a hard convergence faster than DIIS alone. It is the converger used in this tutorial.
  • TRAH (trust-region augmented-Hessian SCF) is the heavy artillery: a genuine second-order method wrapped in a trust region so the Newton step is only taken as far as it can be trusted. It converges cases where DIIS and SOSCF both fail — strongly correlated references, tricky open shells — at a higher per-iteration cost.

The other half of the battle is where you start. A better initial guess puts you inside the basin of the right solution, which cuts iterations and avoids converging to an excited or symmetry-broken state. OpenQP's default is the extended-Hückel guess (type=huckel); it also offers hcore, sap, minao, modhuckel, and JSON-restart guesses. For the full contract behind each option, see the [scf] and [guess] keyword pages.

Input-file style

The runnable deck is inputs/h2o_scf_soscf.inp — water in the 6-31G basis, a spin-restricted PBE reference, converged with the SOSCF converger starting from a Hückel guess. Annotated:

[input]
system=
 8   0.000000000   0.000000000  -0.041061554   # O   (Angstrom)
 1  -0.533194329   0.533194329  -0.614469223   # H
 1   0.533194329  -0.533194329  -0.614469223   # H

charge=0
method=hf              # SCF engine (HF machinery; DFT is switched on by functional)
basis=6-31g
runtype=energy         # single-point SCF energy
functional=pbe         # non-empty -> Kohn-Sham DFT with the PBE functional
d4=false               # no DFT-D4 dispersion correction

[guess]
type=huckel            # extended-Hückel initial orbital guess (the default)
save_mol=true          # keep the converged orbitals/molecule data for restart

[scf]
type=rhf               # closed-shell restricted reference (rhf | uhf | rohf)
multiplicity=1         # singlet ground state
converger_type=soscf   # second-order SCF converger (vs. the default diis / or trah)
maxit=60               # max SCF iterations before giving up
conv=1.0e-6            # energy/density convergence threshold (Hartree)

[dftgrid]
rad_npts=96            # radial quadrature points for the DFT XC grid
ang_npts=302           # angular (Lebedev) points per radial shell
pruned=                # pruning scheme; empty = use the default (unpruned) grid

Key points:

  • The reference is DFT, not bare HF. method=hf selects the SCF machinery, but a non-empty functional (pbe) promotes it to Kohn-Sham DFT. Leaving functional empty would give a Hartree-Fock reference instead.
  • [guess] type chooses where the SCF starts. huckel is the robust default; switching to hcore, sap, or minao changes only the starting orbitals, not the final answer, but can change how many iterations (and whether) you converge. save_mol=true writes the converged orbitals so a later run can restart from them.
  • [scf] converger_type chooses how it gets there. This deck uses soscf. Change one line to switch strategy:
[scf]
converger_type=diis     # default first-order DIIS extrapolation

[scf]
converger_type=trah     # trust-region augmented-Hessian: for hard, stalling cases

All three converge this easy water case to the same energy; the difference shows up on difficult references, where diis may stall and trah still converges. - maxit and conv are the stopping rules. Raise maxit (here 60, above the default of 30) for slow cases; tighten conv (e.g. 1.0e-8) when you need high-precision energies or gradients. - [dftgrid] sets the XC integration grid. rad_npts × ang_npts is a fine (96, 302) grid; a coarse grid can itself cause the SCF energy to jitter and fail to converge, so grid quality and SCF convergence are linked.

Python style

The equivalent calculation with the OpenQP Python API is inputs/h2o_scf_soscf.py. job.theory.dft(...) sets the [input] functional, the basis, and the [scf] reference in one call; any extra keyword (here converger_type, maxit, conv) is forwarded straight into the [scf] section, so it lands exactly where the input-file deck puts it.

from oqp.openqp import OpenQP

job = OpenQP("h2o_scf_soscf", silent=1)

job.molecule(
    """
O   0.000000000   0.000000000  -0.041061554
H  -0.533194329   0.533194329  -0.614469223
H   0.533194329  -0.533194329  -0.614469223
""",
    charge=0,
    multiplicity=1,
)

# DFT/PBE reference; extra kwargs are written to [scf].
# converger_type="soscf" == [scf] converger_type=soscf
job.theory.dft(
    functional="pbe",
    basis="6-31g",
    reference="rhf",         # [scf] type=rhf
    converger_type="soscf",  # second-order SCF instead of DIIS
    maxit=60,                # [scf] maxit
    conv=1.0e-6,             # [scf] conv
)

mol = job.run()

print("SCF energy:", mol.get_scf_energy())
print("Results:", mol.get_results())

To try a different converger, change only the forwarded keyword:

job.theory.dft(functional="pbe", basis="6-31g", reference="rhf", converger_type="diis")   # DIIS
job.theory.dft(functional="pbe", basis="6-31g", reference="rhf", converger_type="trah")   # TRAH

The full script is inputs/h2o_scf_soscf.py.

Run it

Input-file style (from the inputs/ folder):

cd scf-convergence/inputs
openqp h2o_scf_soscf.inp

Python style:

cd scf-convergence/inputs
python h2o_scf_soscf.py

Both need OpenQP installed (pip install openqp) and produce the same SCF energy.

Reading the output

This is a single-point SCF run, so the number you want is the converged SCF (Kohn-Sham) total energy.

  • In the log file (<project>.log) watch the iteration table: each row prints the energy and the density/gradient change, and you want to see that change shrink monotonically past conv=1.0e-6 and the run report converged within maxit iterations. If it hits maxit without converging, that is the signal to switch converger_type (try trah) or improve the [guess].
  • From Python, mol.get_scf_energy() returns the converged SCF total energy, and mol.get_results() returns the results dictionary (the energy field matches what is written to <project>.json).
  • Because [guess] save_mol=true, the converged orbitals are also saved, so a follow-up calculation can restart from them (via a json/auto guess) instead of re-converging from Hückel.

The physical energy is identical no matter which converger you pick — DIIS, SOSCF, and TRAH only differ in how reliably and how fast they reach it. The point of this tutorial is to know which knob to turn when the default DIIS run refuses to converge.

Manual