Self-Consistent Iterative Optimization

• Self-Consistent Iterative Optimization is a method that iteratively updates a state until the induced operator and state become mutually consistent.
• It employs techniques such as damped mixing, DIIS, and quasi-Newton accelerators, along with safeguards like trust-region steps to enhance stability.
• The approach is widely applied in quantum chemistry and electronic structure methods, addressing convergence challenges across various self-consistent frameworks.

Searching arXiv for recent and foundational papers on self-consistent iterative optimization and closely related SCF/Dyson frameworks. Searching "self-consistent iterative optimization SCF Dyson DIIS direct minimization". Self-consistent iterative optimization denotes a class of numerical procedures in which a nonlinear stationarity condition is replaced by a sequence of auxiliary updates, each defined from the current iterate, until a fixed-point, commutator, gradient, or moment-matching residual vanishes. In the sources surveyed here, this pattern appears in finite-temperature Dyson solvers, Kohn–Sham and Hartree–Fock self-consistent field procedures, direct minimization on the Stiefel or Grassmann manifolds, multiconfigurational orbital optimization, tensor-network renormalization, embedding methods based on sampled reduced density matrices, and convex inner–outer schemes for Schrödinger–Poisson systems (Pokhilko et al., 2021, Liu et al., 2013, Cancès et al., 2020, Lacerda-Santos et al., 3 Jul 2025). Although the specific state variable differs—Green’s functions, densities, projectors, orbitals, potentials, or reduced density matrices—the common structure is an iterative enforcement of mutual consistency between a state and the operator, energy model, or environment that the same state induces.

1. Fixed-point and variational structure

A canonical self-consistent formulation is the Kohn–Sham map $K$ acting on the input density $\rho^{\mathrm{in}}$, with self-consistency defined by $K[\rho^\*]=\rho^\*$; in practice one iterates $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$ and terminates when either $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$ or the total-energy change is below a threshold, where $R[\rho]\equiv K[\rho]-\rho$ (Woods et al., 2019). In finite-temperature Green’s-function theory, the Dyson equation is written in imaginary-frequency form as

$$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$

or equivalently

$$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$

so that $G=\mathcal D[G]$ is a nonlinear fixed-point problem because $\Sigma$ is itself a functional of $\rho^{\mathrm{in}}$0 (Pokhilko et al., 2021).

A closely related variational formulation minimizes a smooth functional over an orthogonality-constrained manifold. In one subspace formulation, one seeks

$\rho^{\mathrm{in}}$1

with first-order condition $\rho^{\mathrm{in}}$2, where $\rho^{\mathrm{in}}$3 (Cancès et al., 2020). In direct Kohn–Sham minimization on the Stiefel manifold, one minimizes $\rho^{\mathrm{in}}$4 subject to $\rho^{\mathrm{in}}$5, and the Euler–Lagrange condition is the nonlinear eigenproblem $\rho^{\mathrm{in}}$6 (Zhang et al., 2013). Multiconfigurational formulations retain the same pattern: iCISCF alternates between selected-CI optimization in a current active orbital basis and orbital optimization for the current CI state until both CI and orbital gradients satisfy their thresholds (Guo et al., 2021).

The same abstraction extends beyond orbital theories. In the Non-Linear Helmholtz reformulation of Schrödinger–Poisson, the continuous quantum-electrostatic problem is mapped to

$\rho^{\mathrm{in}}$7

and the exact Schrödinger–Poisson solution is then recovered by an outer loop that updates the integrated local density of states and re-solves the inner convex problem (Lacerda-Santos et al., 3 Jul 2025). In reduced-density-matrix embedding, one repeatedly constructs a correlated subspace Hamiltonian from environment parameters, samples one- and two-body reduced density matrices on the subspace, and updates the classical environment so that energy-weighted moments or orbital gradients become self-consistent (Tilly et al., 2021).

Problem class	Self-consistency variable	Auxiliary update
Kohn–Sham SCF	$\rho^{\mathrm{in}}$8 or $\rho^{\mathrm{in}}$9	Linear eigenproblem for $K[\rho^\*]=\rho^\*$0 (Liu et al., 2013)
Finite-$K[\rho^\*]=\rho^\*$1 Dyson	$K[\rho^\*]=\rho^\*$2	Dyson update $K[\rho^\*]=\rho^\*$3 (Pokhilko et al., 2021)
Projector minimization	$K[\rho^\*]=\rho^\*$4	Retraction after tangent-space step (Cancès et al., 2020)
Schrödinger–Poisson via NLH	$K[\rho^\*]=\rho^\*$5	Convex minimization of $K[\rho^\*]=\rho^\*$6 (Lacerda-Santos et al., 3 Jul 2025)
CASSCF / iCISCF	Orbitals and CI coefficients	Alternating CI and orbital optimization (Guo et al., 2021)

Taken together, these formulations show that “self-consistency” is not restricted to Roothaan-type orbital updates. A plausible implication is that the defining feature is not the linear eigenproblem itself, but the repeated enforcement of compatibility between a state and the effective operator or environment generated by that state.

2. Residuals, objective functions, and stationarity measures

The residual chosen to monitor self-consistency determines both diagnostics and acceleration. In Kohn–Sham SCF, the standard residual map is $K[\rho^\*]=\rho^\*$7 (Woods et al., 2019). In the finite-temperature Dyson problem, two residuals were compared: the difference residual

$K[\rho^\*]=\rho^\*$8

and the commutator residual

$K[\rho^\*]=\rho^\*$9

Because the Dyson equation is exactly $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$0, the commutator vanishes if and only if $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$1 solves it, whereas $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$2 only tracks the iteration history; empirically, the commutator residual was more sensitive to violations of the Dyson equation and yielded dramatically improved convergence behavior (Pokhilko et al., 2021).

Commutator residuals also appear in conventional SCF. In the iterative stochastic subspace SCF method, convergence is monitored by the Frobenius norm

$$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$3

together with the energy change $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$4 (Loos et al., 2013). OpenOrbitalOptimizer uses the commutator error $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$5, which is zero at convergence, as the DIIS error vector (Lehtola et al., 29 Mar 2025). In projector optimization, the Riemannian gradient is

$$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$6

so stationarity again reduces to a commutator-like tangential residual (Cancès et al., 2020).

When the method is recast as explicit minimization, the relevant quantity is an objective functional rather than a residual norm. For the NLH equation, the convex functional

$$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$7

has gradient $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$8 and Hessian $$\rho^{\mathrm{in}}_{n+1}=f(\{\rho^{\mathrm{in}}_i,\rho^{\mathrm{out}}_i\}_{i\le n})$$9; since $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$0 is symmetric positive semi-definite and $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$1, $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$2 is globally convex (Lacerda-Santos et al., 3 Jul 2025). In orthogonal CCA, the subproblem is written as a trace-fractional maximization on the Stiefel manifold, and the customized SCF iteration is built from the KKT system of that objective (Zhang et al., 2019).

The same structural distinction appears in SCF line search. Adaptive damping evaluates the dual energy functional $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$3 and accepts a step if either the trial energy decreases or the preconditioned residual norm decreases; an Armijo criterion is also available (Herbst et al., 2021). This suggests that self-consistent iterative optimization admits two complementary views: a fixed-point view based on a residual map, and a variational view based on energy decrease or manifold stationarity. The strongest algorithms in the surveyed literature often combine both.

3. Update rules, subspace accelerators, and safeguarded steps

The simplest self-consistent update is damped mixing. In Kohn–Sham DFT this takes the form

$\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$4

or, in potential form,

$\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$5

where $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$6 is a preconditioner (Woods et al., 2019, Herbst et al., 2021). Fixed damping is robust under suitable assumptions but often slow, so most practical solvers enlarge the update space with multisecant or quasi-Newton information.

Pulay’s DIIS is the standard subspace construction. One stores recent iterates and residuals and chooses coefficients $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$7 that minimize the norm of the extrapolated residual subject to $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$8 (Woods et al., 2019). In Dyson solvers, the extrapolated object may be $\|R[\rho^{\mathrm{in}}_n]\|\ll\mathrm{tol}$9, with residual overlaps computed in the time domain via trace integrals over $R[\rho]\equiv K[\rho]-\rho$0; the corresponding linear system is the standard constrained least-squares DIIS system (Pokhilko et al., 2021). OpenOrbitalOptimizer implements DIIS together with EDIIS, ADIIS, and the optimal damping algorithm; EDIIS minimizes the Hartree–Fock energy evaluated at an interpolated density under simplex constraints, while ADIIS uses a quadratic expansion around the current density (Lehtola et al., 29 Mar 2025).

Beyond DIIS, several distinct accelerators appear. KAIN forms an inexact Newton step in the Krylov subspace spanned by differences of past iterates and residuals, though in the Dyson setting it must use a difference-type residual because the commutator residual is orthogonal to the state vector and leads to a singular system (Pokhilko et al., 2021). Broyden and multisecant Broyden updates approximate either the Jacobian or inverse Jacobian of the SCF map, and the SCF-$R[\rho]\equiv K[\rho]-\rho$1 benchmark found Broyden II + Kerker and Pulay DIIS + Kerker to be among the strongest standard methods (Woods et al., 2019). QUOTR replaces the Roothaan–Hall step by unitary orbital rotations generated by a preconditioned L-BFGS algorithm with trust-region restriction, thereby enforcing orthonormality through $R[\rho]\equiv K[\rho]-\rho$2 with $R[\rho]\equiv K[\rho]-\rho$3 (Slattery et al., 2023).

Safeguards are pervasive. In Dyson solvers, step restriction rescales the update coefficients so that $R[\rho]\equiv K[\rho]-\rho$4, with $R[\rho]\equiv K[\rho]-\rho$5–3 used as a trust radius (Pokhilko et al., 2021). In adaptive damping for SCF, a backtracking line search adjusts $R[\rho]\equiv K[\rho]-\rho$6 using an inexpensive quadratic model of $R[\rho]\equiv K[\rho]-\rho$7, and the model-quality ratio determines whether to accept the model minimizer or simply halve the damping (Herbst et al., 2021). Heavy-ball dynamics augments the orbital update with a momentum term $R[\rho]\equiv K[\rho]-\rho$8, while Laplacian-based potential preconditioning applies $R[\rho]\equiv K[\rho]-\rho$9 to suppress high-$$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$0 oscillations in the mean-field potentials (Ryssens et al., 2018).

A recurrent theme is that efficient self-consistent optimization usually requires both acceleration and stabilization. The literature does not present a single dominant mechanism; rather, it shows repeated combinations of extrapolation, quasi-Newton information, trust regions, line search, and preconditioning.

4. Convergence theory, spectral gaps, and failure modes

The strongest available theory is local and assumption-dependent. For Kohn–Sham SCF, one analysis proves global convergence from an arbitrary initial point and local linear convergence from an initial point sufficiently close to the solution under two stringent assumptions: a fixed spectral gap $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$1 between occupied and unoccupied states at every iterate, and a uniform upper bound $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$2 on the second derivatives of the exchange–correlation functional (Liu et al., 2013). Under these assumptions, the basic SCF iteration decreases the total energy monotonically and converges; locally, the contraction factor is

$$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$3

(Liu et al., 2013). The same paper interprets SCF as an inexact Newton method that replaces the full Hessian $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$4 by the Hamiltonian block $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$5, with the neglected term $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$6 governing convergence behavior (Liu et al., 2013).

A complementary analysis on the manifold of rank-$$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$7 orthogonal projectors compares damped SCF and Riemannian gradient descent. Under $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$8 smoothness and nondegeneracy, gradient descent converges linearly with asymptotic rate $$G^{-1}(i\omega_n)=G_0^{-1}(i\omega_n)-\Sigma[G](i\omega_n),$$9, while damped SCF converges linearly under a strict gap $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$0 with rate $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$1 (Cancès et al., 2020). The comparison is explicit: SCF is highly sensitive to the gap $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$2, whereas gradient descent is sensitive to the overall spectral span and not directly to the small gap (Cancès et al., 2020).

Globalized variants can relax the practical fragility of plain SCF. Backtracking Armijo damping guarantees monotonic energy decrease, residual convergence to zero, and global convergence of $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$3 under the paper’s assumptions; once close enough to a solution with positive-definite Hessian, local linear or superlinear convergence is recovered depending on Anderson acceleration (Herbst et al., 2021). In stochastic SCF, convergence in the mean-square sense or in probability can be proved when the underlying fixed-point map is locally contractive and the stochastic error satisfies boundedness or second-moment assumptions (Ko et al., 2021).

The surveyed literature also gives explicit examples of failure. Charge sloshing in metals amplifies long-wavelength density errors because the Jacobian of the Kohn–Sham map has small-$$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$4 eigenvalues of order $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$5 (Woods et al., 2019). Bare SCF can stall or oscillate on small-gap or metallic systems (Zhang et al., 2013). In a nonquadratic subspace optimization model, the simple undamped SCF map undergoes a period-doubling cascade and enters a chaotic regime beyond a critical parameter $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$6, demonstrating that SCF can fail by genuine chaos rather than simple oscillation (Cancès et al., 2020). In the p-Laplacian eigenproblem, convergence slows as $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$7 because the eigenvalue gap around the target eigenvalue shrinks, and regularization by “softabs” is needed to avoid ill-conditioning (Upadhyaya et al., 2021).

These results establish a consistent picture: self-consistent iterative optimization is efficient when the effective map is contractive or the neglected curvature is controlled, but it can become arbitrarily slow or unstable in small-gap, strongly nonlinear, or ill-conditioned regimes.

5. Multilevel, embedding, and nonstandard realizations

Several methods depart from classical SCF while retaining the same self-consistent logic. Direct minimization on the Stiefel manifold updates

$$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$8

using QR or WY retractions and often a Barzilai–Borwein initialization followed by a non-monotone Armijo line search (Zhang et al., 2013). This avoids repeated large eigenvalue problems and inherits global convergence from standard manifold-optimization theory (Zhang et al., 2013). QUOTR pushes this further with unitary rotations, L-BFGS inverse-Hessian approximation, and a low-rank trust-region solve (Slattery et al., 2023).

Multilevel decomposition also appears in fragment and active-space methods. The iOI approach uses an outer macroiteration over fragment mergers and an inner microiteration of SCF solves on each subsystem, with subsystem growth governed by a tail-population criterion; the final global step uses block-diagonalization with orbital freezing (Wang et al., 2021). iCISCF alternates between configuration selection and orbital optimization, using Jacobi rotation for active–active updates and quasi-Newton BFGS updates for the remaining orbital rotations; after variational convergence it applies an Epstein–Nesbet PT2 correction on the residual space (Guo et al., 2021). In both cases, self-consistency is distributed across coupled subproblems rather than concentrated in a single density-mixing loop.

Tensor-network and embedding formulations exhibit the same pattern. In finite-temperature PEPS, every Suzuki–Trotter step increases the bond dimension, after which one truncates back to the original $$G(i\omega_n)=G_0(i\omega_n)+G_0(i\omega_n)\,\Sigma[G](i\omega_n)\,G(i\omega_n),$$9, computes the full environment of the truncated tensor network, updates the bond isometry, and repeats until self-consistency $G=\mathcal D[G]$0 is reached (Czarnik et al., 2014). The update of the isometry maximizes an approximate fidelity functional, and because each update is optimal at fixed environment, $G=\mathcal D[G]$1 never decreases (Czarnik et al., 2014). In reduced-density-matrix-sampling embedding, the outer loop constructs a subspace Hamiltonian, solves it on a quantum device, samples one- and two-body reduced density matrices, and updates orbital rotations or auxiliary parameters so that gradients or energy-weighted moments match between the correlated subspace and the classical environment (Tilly et al., 2021).

The NLH scheme for Schrödinger–Poisson makes the multilevel structure especially explicit. An inner convex minimization solves the approximate NLH problem, while an outer loop updates the true local density of states and re-solves the convex problem until the exact self-consistent Schrödinger–Poisson solution is recovered (Lacerda-Santos et al., 3 Jul 2025). The paper states that empirically only one or two outer updates are required to recover the exact self-consistent solution to machine precision (Lacerda-Santos et al., 3 Jul 2025). This suggests that “self-consistent iterative optimization” can encompass nested architectures in which inner problems are made provably well behaved and outer loops restore the original nonlinear physics.

6. Empirical behavior, benchmark results, and practical prescriptions

The empirical literature emphasizes both performance variability and transferable heuristics. In finite-temperature Dyson solvers, commutator-based CDIIS and LCIIS outperform difference-residual schemes for all considered molecular and solid systems within both GW and GF2 (Pokhilko et al., 2021). For the beryllium atom in GW at $G=\mathcal D[G]$2, plain damping with $G=\mathcal D[G]$3 converges at low temperature but diverges at $G=\mathcal D[G]$4, DDIIS often oscillates or diverges at high temperature, KAIN behaves similarly to DDIIS, whereas CDIIS and LCIIS converge in 6–10 iterations for all tested $G=\mathcal D[G]$5 (Pokhilko et al., 2021). For solids such as Si, NiO, and BiVO$G=\mathcal D[G]$6, CDIIS and LCIIS reach $G=\mathcal D[G]$7 and $G=\mathcal D[G]$8 in 8–15 iterations, with wall-clock speedup $G=\mathcal D[G]$9 relative to damping, and extra overhead $\Sigma$0 of a single GW/GF2 iteration (Pokhilko et al., 2021).

In Kohn–Sham DFT, the SCF-$\Sigma$1 benchmark suite of about 56 inputs provides robustness $\Sigma$2 and efficiency $\Sigma$3 metrics for algorithm comparison (Woods et al., 2019). Unpreconditioned linear mixing with $\Sigma$4 converged about $\Sigma$5 of cases with $\Sigma$6; Pulay DIIS + Kerker with $\Sigma$7, $\Sigma$8, and $\Sigma$9 reached $\rho^{\mathrm{in}}$00 at $\rho^{\mathrm{in}}$01; “Periodic Pulay” and “Restarted Pulay” improved robustness to about $\rho^{\mathrm{in}}$02–$\rho^{\mathrm{in}}$03 with only about $\rho^{\mathrm{in}}$04 loss in efficiency; ensemble-DFT converged about $\rho^{\mathrm{in}}$05 of cases but with $\rho^{\mathrm{in}}$06 (Woods et al., 2019). The recommended default in that study is preconditioned multisecant Pulay/DIIS with history $\rho^{\mathrm{in}}$07, $\rho^{\mathrm{in}}$08, and $\rho^{\mathrm{in}}$09, together with Kerker preconditioning and occasional history refreshes (Woods et al., 2019).

Adaptive damping was designed to remove manual tuning of $\rho^{\mathrm{in}}$10. On elongated Al supercells, adaptive damping with default $\rho^{\mathrm{in}}$11 converged as quickly as the best fixed $\rho^{\mathrm{in}}$12, typically in 11–17 diagonalizations for Al$\rho^{\mathrm{in}}$13 and 12–21 for Al$\rho^{\mathrm{in}}$14; on GaAs it automatically down-scaled to avoid large error spikes caused by nonlinear Anderson artifacts; on transition-metal systems it converged all but one case, and increasing the minimum damping to $\rho^{\mathrm{in}}$15 recovered the last (Herbst et al., 2021). In well-behaved regimes, the step is accepted at the first trial, so the method uses only one diagonalization per SCF iteration, exactly as in fixed damping; in the worst case, backtracking can introduce a small constant-factor overhead (Herbst et al., 2021).

Alternative orbital optimizers exhibit different trade-offs. QUOTR required median 16–20 Fock builds on standard G2 test sets, compared with 12–15 for DIIS, but had zero local-minimum failures whereas DIIS had 3–8 and exact Newton or augmented-Hessian methods had 1–4 (Slattery et al., 2023). For the 1PLW peptide with vanishing Kohn–Sham gap, DIIS failed for LDA and PBE whereas QUOTR converged all functionals and found gaps as small as $\rho^{\mathrm{in}}$16 (Slattery et al., 2023). Heavy-ball dynamics plus potential preconditioning reduced total SCF iterations by a factor of $\rho^{\mathrm{in}}$17–$\rho^{\mathrm{in}}$18 in Skyrme-EDF calculations and required essentially no manual fine-tuning beyond isolated parametrizations such as UNEDF1 (Ryssens et al., 2018).

Several practical prescriptions recur across otherwise distinct domains. Always preferring commutator residuals over simple difference residuals is recommended for finite-$\rho^{\mathrm{in}}$19 Dyson solvers (Pokhilko et al., 2021). Small subspaces, often dimension $\rho^{\mathrm{in}}$20–$\rho^{\mathrm{in}}$21, are favored for robustness in CDIIS and LCIIS (Pokhilko et al., 2021). Trust-radius or trust-region control is repeatedly invoked to prevent overshooting (Pokhilko et al., 2021, Slattery et al., 2023). Warm starts, cooling schemes, p-continuation, finite-$\rho^{\mathrm{in}}$22 smearing, and ensemble formulations are recommended for pathological cases with stretched bonds, multiple open shells, or low spectral gaps (Pokhilko et al., 2021, Upadhyaya et al., 2021, Woods et al., 2019). The common lesson is that the most reliable self-consistent optimizers do not rely on a single mechanism; they combine a physically informative residual or objective, a subspace or quasi-Newton accelerator, and an explicit safeguard against noncontractive steps.