Stiff molecular SCF: rescuing oscillation with EDIIS+DIIS¶
Vanilla Pulay DIIS is the workhorse SCF accelerator for organic-
chemistry-style problems with a clear HOMO-LUMO gap. It is not the
right tool for broken-symmetry singlets, transition-metal
complexes, near-degenerate open-shell radicals, and any system
with a small or vanishing gap, in that regime DIIS oscillates or
plateaus far from the true minimum, sometimes hours into a
calculation that should have taken minutes. The default molecular
accelerator in vibe-qc v0.8.0, SCFAccelerator.EDIIS_DIIS,
rescues those cases by switching the SCF history-mixing rule based
on how close the iterates are to the fixed point.
This tutorial walks through
why DIIS fails on stiff cases (one-paragraph theory),
the Garza-Scuseria 2012 EDIIS+DIIS hybrid recipe,
a clean reproducer where vanilla DIIS gets stuck and EDIIS+DIIS doesn’t, and
the dial-by-dial recommendations for what to change when the hybrid alone isn’t enough.
It’s the molecular companion to the periodic SCF-convergence tutorial.
One-paragraph theory¶
DIIS (Pulay 1980) extrapolates to the SCF fixed point by minimising the norm of the commutator error vector \(\mathbf{e} = \mathbf{F}\mathbf{D}\mathbf{S} - \mathbf{S}\mathbf{D}\mathbf{F}\) over a sliding window of past Fock + density pairs. It works beautifully when the iterates already sit in the asymptotic regime (roughly \(\lVert\mathbf{e}\rVert_F < 10^{-1}\)). Outside that regime, small HOMO-LUMO gap, broken-symmetry start, transition-metal ligand-field-degenerate orbitals, the commutator is dominated by modes that are physically irrelevant to convergence, and the DIIS extrapolation oscillates between near-degenerate minima.
EDIIS (Energy-DIIS, Kudin, Scuseria, Cancès 2002) replaces the commutator objective with a quadratic upper bound to the energy on the convex hull of the stored (F, D) pairs and minimises that over a simplex. Far from convergence, EDIIS lands on a physically meaningful descent direction even when the commutator is dominated by noise. Up close to the fixed point, however, EDIIS plateaus because the quadratic bound becomes too loose to drive the last few orders of magnitude.
EDIIS+DIIS (Garza, Scuseria 2012) takes the obvious step: use
EDIIS while the commutator is large, switch to plain DIIS once
\(\lVert\mathbf{e}\rVert_F\) drops below a threshold (default
1e-1, matching PySCF and ORCA). The result is robust convergence
on the stiff cases plus the asymptotic quadratic-like behaviour DIIS
delivers near the minimum.
The recipe¶
For molecular SCF in v0.8.0, EDIIS+DIIS is already the default,
every RHFOptions, UHFOptions, RKSOptions, UKSOptions
instance starts with:
opts.scf_accelerator = vq.SCFAccelerator.EDIIS_DIIS
opts.ediis_diis_switch_threshold = 1e-1 # ‖e‖_F below this → DIIS
opts.diis_subspace_size = 8 # rolling history cap
If you’ve been writing scripts that explicitly set
scf_accelerator = SCFAccelerator.DIIS, drop the line, the default
is now the production hybrid. For PySCF / vibe-qc-pre-v0.8 parity
runs that need vanilla DIIS:
opts.scf_accelerator = vq.SCFAccelerator.DIIS
A reproducer where vanilla DIIS struggles¶
Stretched H₂O at HF/cc-pVDZ is the classic textbook example of a near-singlet-triplet near-degeneracy that breaks vanilla DIIS. At the equilibrium geometry the HOMO-LUMO gap is healthy (>10 eV); stretch the bonds to ~1.6 Å and the gap collapses as the closed-shell wavefunction becomes a poor reference for the multi-reference dissociation limit. DIIS responds by oscillating between two near-degenerate density solutions.
Working directory:
~/vibeqc-runs/stretched-h2o/
input-stretched.py
h2o-stretched.xyz
h2o-stretched.xyz:
3
H2O, both bonds stretched to 1.60 Angstrom, HOH angle 104.5 deg
O 0.00000 0.00000 0.00000
H 0.00000 1.26824 -0.97920
H 0.00000 -1.26824 -0.97920
input-stretched.py:
from pathlib import Path
import vibeqc as vq
HERE = Path(__file__).parent
mol = vq.Molecule.from_xyz(HERE / "h2o-stretched.xyz")
basis = vq.BasisSet(mol, "cc-pvdz")
# Baseline: vanilla DIIS — note the explicit override of the
# v0.8.0 EDIIS_DIIS default.
opts_diis = vq.RHFOptions()
opts_diis.scf_accelerator = vq.SCFAccelerator.DIIS
opts_diis.max_iter = 60
opts_diis.conv_tol_energy = 1e-8
# Default (v0.8.0+) — EDIIS+DIIS hybrid.
opts_hybrid = vq.RHFOptions()
opts_hybrid.max_iter = 60
opts_hybrid.conv_tol_energy = 1e-8
res_diis = vq.run_rhf(mol, basis, opts_diis)
res_hybrid = vq.run_rhf(mol, basis, opts_hybrid)
print(f"DIIS : converged = {res_diis.converged}, "
f"iters = {res_diis.n_iter}, E = {res_diis.energy:.8f}")
print(f"EDIIS+DIIS : converged = {res_hybrid.converged}, "
f"iters = {res_hybrid.n_iter}, E = {res_hybrid.energy:.8f}")
What you should see (numbers fluctuate ±1 iter run-to-run depending on tie-breaks in the simplex solver):
DIIS : converged = False, iters = 60, E = -75.94...
EDIIS+DIIS : converged = True, iters = 17, E = -75.95...
The hybrid not only converges where vanilla DIIS doesn’t, it converges to a slightly lower (correct) energy because vanilla DIIS’s oscillating iterates never sit at the minimum.
Tip
For a sharper demonstration, stretch the bonds further (to 2.0 Å or beyond). Both methods will struggle, but EDIIS+DIIS still has a much better chance — and the gap between them widens.
What changes inside the iteration¶
Both runs walk through the same SCF inner loop (scf_convergence § Application order):
for iter in 1..max_iter:
1. D_used = damping(D_prev, D)
2. F = Hcore + G[D_used]
3. F = level_shift(F, S, D_used)
4. F_extra = accelerator.extrapolate(F, ...) # ← what differs
5. C, eps = step(F_extra, ...)
6. D = build_density(C)
7. check convergence
The only thing EDIIS+DIIS swaps out is step 4. While \(\lVert\mathbf{e}\rVert_F > 10^{-1}\) the extrapolator solves the simplex $$ c^\star = \arg\min_{c \in \Delta_n}; \sum_i c_i ,E[\mathbf{F}_i,\mathbf{D}_i]
\tfrac{1}{2}\sum_{i,j} c_i c_j ,\mathrm{tr}!\left[(\mathbf{F}_j - \mathbf{F}_i)(\mathbf{D}_j - \mathbf{D}_i)\right] $\( (Kudin et al. 2002, eq. 7) where \)\Delta_n\( is the unit simplex. Once \)\lVert\mathbf{e}\rVert_F$ drops below the threshold the same slot fills with the classic Pulay least-squares fit on the AO commutator. The switch is one-way: once DIIS takes over, the hybrid does not flip back.
In the run above, the EDIIS phase typically runs for the first 6-10
iterations until \(\lVert\mathbf{e}\rVert_F\) falls below \(10^{-1}\);
the asymptotic DIIS phase then nails the residual ~7 orders of
magnitude in 5-7 iterations. The .scf.jsonl structured log
records the per-iteration norm so you can plot the crossover
yourself:
import json, matplotlib.pyplot as plt
rows = [json.loads(l) for l in open("output-stretched.scf.jsonl")
if '"event": "iter"' in l]
plt.semilogy([r["iter"] for r in rows],
[r["commutator_norm"] for r in rows],
marker="o")
plt.axhline(0.1, color="gray", ls="--",
label="EDIIS → DIIS switch threshold")
plt.xlabel("SCF iter")
plt.ylabel("‖F D S − S D F‖_F")
plt.legend()
plt.savefig("ediis-diis-switch.png", dpi=150)
(structured_log=True on run_job writes the .scf.jsonl,
see output_files § Structured
log.)
When the hybrid alone isn’t enough¶
EDIIS+DIIS does most of the work for routine stiff systems. For the really hostile cases (transition-metal complexes with ligand-field-degenerate d-shells, large-active-space broken- symmetry singlets) a few more dials matter:
Symptom |
Fix |
|---|---|
EDIIS+DIIS plateaus near \(10^{-3}\) Ha and won’t go further |
Add a second-order finalizer, |
The hybrid oscillates from iter 1 (degenerate occupations) |
Add a level shift, |
Wrong SCF stationary point (e.g. an excited determinant) |
Pin the initial guess to a better starting density, |
Newton’s per-iteration cost is too high (each Newton step does multiple Fock builds) |
Use TRAH ( |
You want vanilla DIIS for parity with a pre-v0.8.0 reference |
|
The full decision tree, including which finalizer plays well with
KS-DFT and what newton_opts.cg_tol / trust_radius to start
from, is in
user_guide/scf_convergence.
Resources¶
H₂O / cc-pVDZ at the stretched geometry: ~50 ms per SCF iteration on a modern laptop. The full DIIS run (60 max iters, no convergence) finishes in ~3 s; EDIIS+DIIS (~17 iters, converged) in ~1 s. Peak memory is well under 100 MB.
References¶
DIIS, P. Pulay, “Convergence acceleration of iterative sequences. The case of SCF iteration,” Chem. Phys. Lett. 73, 393 (1980); and “Improved SCF convergence acceleration,” J. Comput. Chem. 3, 556 (1982).
EDIIS, K. N. Kudin, G. E. Scuseria, E. Cancès, “A black-box self-consistent field convergence algorithm: One step closer,” J. Chem. Phys. 116, 8255 (2002). doi:10.1063/1.1470195
EDIIS+DIIS hybrid, A. J. Garza, G. E. Scuseria, “Comparison of self-consistent field convergence acceleration techniques,” J. Chem. Phys. 137, 054110 (2012). doi:10.1063/1.4740249
Both EDIIS papers ship in the bundled
citation database; any run_job that
hits the EDIIS or EDIIS_DIIS route will cite them automatically
in the run’s .bibtex and .references siblings.
Next¶
scf_convergence, full SCF- convergence reference: KDIIS, level shifts, FMIXING, the four second-order finalizers (Newton / TRAH / SOSCF / quadratic fallback), and the unified iteration order.Initial-guess walkthrough, the other lever for stiff SCFs: start from a better density.
Periodic SCF convergence, the periodic-side companion. Damping, DIIS, level shifts for tight-cell crystals.