The SCF loop every Hartree-Fock method runs

Hartree-Fock replaces the electron-electron repulsion by the average field of the other electrons. That field depends on the orbitals, and the orbitals depend on the field, so the equations are solved by iteration: guess a density, build the Fock matrix from it, diagonalize, occupy the lowest orbitals, and repeat until the density stops changing. The converged density is a fixed point of that map, and a fixed point is also a stationary point of the energy: F D S - S D F = 0.

scf_iterate is that loop, written once for RHF (one total density), UHF (a density per spin) and ROHF (a density per spin, one set of orbitals). One iteration, with X^T S X = I (m columns kept):

1. J, K from the densities; the Fock matrices (F = H + J - K/2 for RHF, F_s = H + J_a + J_b - K_s otherwise).
2. The energy, E = E_nuc + (1/2) sum_d Tr D_d (H + F_d).
3. The commutator in the orthonormal basis: F' P' - P' F' per density for RHF and UHF; for ROHF the effective matrix of lib/hf/coupling and its commutator with P_c' + P_o'/2.
4. Stop when the energy changed less than e_tol and the largest commutator element is below grad_tol. The density returned is the one this iteration started from, so the energy, the density and the error all describe the same state.
5. Otherwise DIIS picks the matrix to diagonalize (from the iteration after diis_start on; before that the new Fock matrix as it is), a level shift (when one is set) raises the virtual block, and the lowest orbitals give the next densities.

Each iteration allocates what it computes and keeps only the new density; a compiled loop frees the rest when the iteration ends, so 200 iterations run in the memory of one [test "200 iterations run in fixed memory and stay on the solution"] (lib/hf/rhf).

The stages each have several members, picked by a field of the run's ScfOptions (lib/hf/opts): the guess (gs_of), the orthogonalizer (orth_of), the eigensolver (eig_of), the accelerator (accel_extrapolate), the ROHF coupling and the J/K builder (jk_of). Each of those calls its members by name, so the default options run the production members directly; the registries that list the members are for the tests and sweeps that run them side by side.

One entry point. scf_solve(setup, start) runs every method from spin densities start[2, n, n]: RHF adds them, UHF and ROHF iterate on them, and ROHF with fewer alpha than beta electrons runs with the counts swapped and swaps the densities back. The three ways in are thin:

scf(setup)                     from the guess setup.options.guess
scf_run_from(setup, Da, Db)    from given spin densities (import mode)
scf_step(setup, Da, Db)        one plain step: F from D, diagonalize, occupy

scf(setup) is what lib/hf/scf's rhf, uhf, rohf and scf_with run after integrals(mol, o) has built the setup.

At a converged state scf_step gives the densities back, which checks the Fock build, the diagonalization and the occupation rule against PySCF one step at a time, whatever path an iteration took.

The result is an ScfResult: the energy, whether it converged and after how many iterations, the densities and orbitals (total, then per spin in alpha and beta), the largest |F D S - S D F| in the AO basis (grad), <S^2> (0 for RHF) and the energy by a second formula (energy_check, E_nuc + Tr D H + (1/2) Tr D J - (1/4) Tr D K for RHF). For RHF the top-level orbitals are the spatial ones with occupations 2 and 0, and each spin holds half the density. For ROHF they are the eigenvectors of the effective matrix with occupations 2, 1 and 0. For UHF they are the alpha orbitals. fock is the Fock matrix for RHF and the average of the spin Fock matrices otherwise.

Used with import hf.scf_loop.