Speeding up the SCF fixed point by DIIS, plain iteration and damping

One SCF iteration maps a Fock matrix to a density (diagonalize, occupy) and the density back to a Fock matrix (J and K). Plain iteration applies that map until nothing changes. It converges when the map is a contraction near the solution, which is often slowly and sometimes not at all: the density can swing between two shapes forever.

An accelerator chooses the Fock matrix the next step diagonalizes. The members:

• accel_cdiis (member 0, production): Pulay's DIIS with the commutator error. At a solution F P - P F = 0 (in the orthonormal basis), so the commutator e_i of each stored step measures how far that step is from it. DIIS looks for the combination F = sum c_i F_i, sum c_i = 1, whose combined error sum c_i e_i is smallest, which is a small linear system:

[ B   -1 ] [ c      ]   [  0 ]
[ -1   0 ] [ lambda ] = [ -1 ],      B_ij = <e_i, e_j>.

If the errors depended linearly on the Fock matrix this would be exact after as many steps as the error has independent directions (it is then GMRES on the linearized map); the test below checks that on a linear problem.

• accel_none: plain iteration, the newest Fock matrix as it is.
• accel_damp: F = (1 - a) F_new + a F_prev, where F_prev is the matrix the previous step used, with a = opts[OPT_DAMP] (ScfOptions.damping, 0.3). Damping cannot speed up a slow contraction, but it can turn an oscillation into a convergent spiral.

State. The history is the record Diis, made once before the loop by accel_new(k, m, opts) for k matrices of m x m a step:

hist[2 nd + 1, k, m, m]   slots 0..nd-1   Fock matrices (ring)
                          slots nd..2nd-1 their errors (same ring)
                          slot 2 nd       the matrix the last step used
meta[ACCEL_META_LEN]      count stored, newest slot, capacity nd, steps
gram[nd, nd]              <e_s, e_t> for ring slots s and t

with nd = opts[OPT_DIIS_N]. accel_push(d, F, E) writes the newest pair into the ring and adds the new error's inner products with the stored ones to gram, k inner products a step instead of the k (k + 1) / 2 a rebuilt B would need (PySCF's lib/diis.py keeps its B the same way). accel_extrapolate(d, opts) runs the member opts[OPT_ACCEL], records its matrix in slot 2 nd and returns it. The ring is the accelerator's state, filled in place; it does not grow.

Member signature: (hist, meta, opts) -> [ndens, m, m]. A member called through the registry gets no gram, so accel_cdiis recomputes it from the history with the same inner product (accel_gram) and gives the same bits; accel_extrapolate calls accel_cdiis_gram with the kept one. With two densities (UHF), DIIS uses one set of coefficients for both spins and adds the two spins' error products in B, as PySCF does.

Used with import hf.accelerate.