06Making the classical solution trustworthy

A successful-looking energy is the beginning of an argument

Our nitrogen model has not changed: neutral singlet N₂ at 1.85 Å, initially in spherical cc-pVDZ. Lesson 5 explained how SCF produces a determinant and energy. Now we ask whether that result is a sufficiently accurate solution of the intended mean-field problem, and whether another nearby determinant would lower its energy.

These are different questions. Nonconvergence means the numerical iteration has not satisfied its stopping conditions. Instability means an apparently stationary solution has a nearby downhill direction in an allowed orbital space. A metastable solution is a local minimum that may not be the lowest minimum accessible elsewhere. Model error remains even for a perfectly converged, stable HF determinant because HF and a finite basis do not describe all electronic correlation. Printing more decimal places addresses none of those distinctions.

Everything below is source-grounded preparation. No SCF calculation, stability check, or test was executed. Small mathematical examples are explicitly illustrative. We will not assign the stability results of the repository’s other nitrogen geometries to our 1.85-Å cc-pVDZ case.

Why the same model can produce difficult iterations

SCF alternates between an electron density and the field that density generates. Near a well-behaved solution, a small density error may shrink each time around the loop. Near a small occupied–virtual energy gap, a small field change can strongly rotate orbitals or reorder occupations. The next density may then undo the previous change. This gives oscillation, slow convergence, or large jumps between electronic solutions.

An illustrative one-variable iteration makes the numerical issue visible. Let $x_k$ represent one component of the error at step $k$, not a nuclear coordinate. If locally $x_{k+1}=-1.2x_k$, the sign flips and the magnitude grows. Mixing half of the old and half of the proposed new value instead gives $x_{k+1}=0.5x_k+0.5(-1.2x_k)=-0.1x_k$, which shrinks. This example explains why damping can help. It does not assert that the native solver exposes a high-level damping keyword or that every molecular SCF problem reduces to this scalar map.

A poor initial guess can land in a different basin of attraction. A symmetry-preserving initial guess can stay on a symmetry-preserving stationary branch even when a symmetry-breaking solution is lower. An optimizer cannot always discover a downhill direction whose first derivative is exactly zero at the symmetric starting point. That is why “use more iterations” and “check stability” are not interchangeable remedies.

Start with a meaningful guess and a consistent problem

The native SCF adapter selects its atom-based initial-density path for a basis name or explicit BasisSet, and its user-provided-density path for Orbitals. The atom-guess code reads available atomic-guess data or generates atomic densities and assembles them on the molecule. Generated guesses use spherically averaged atomic HF. This gives a sensible starting electron distribution before molecular bonding has been optimized.

For Orbitals, the adapter builds a density from occupied coefficient columns using the requested alpha/beta populations. Reusing a previous solution can therefore be much better than starting afresh. It must still describe the same compatible basis and centers. Passing cc-pVDZ coefficients does not automatically project them into cc-pVTZ, and a stale basis object does not move its centers merely because a different Structure is passed beside it.

The low-level engine contains other initial-density strategies, including core-Hamiltonian and SOAD paths. These are not all exposed as settings on create("scf_solver", "qdk"). The supported high-level control here is the input choice. The PySCF adapter similarly uses a compatible orbital input to make an initial density, otherwise calls PySCF’s default kernel() behavior. It does not expose every PySCF initial-guess option through its small QDK settings schema.

For closed-shell unrestricted requests, both adapters warn that automatic spin-symmetry breaking is not supplied. A restricted guess copied into two equal spin channels can simply return the restricted solution inside the larger unrestricted space. A deliberately perturbed unrestricted guess has a different purpose from an arbitrary random starting density: it probes an identified electronic direction while retaining the molecule and electron count.

Read the stopping conditions, not just the last energy change

For ordinary native RHF/UHF, the orbital stationarity residual is associated with

$F$ is the Fock matrix, $P$ the density matrix, and $S$ the AO overlap matrix. A residual measures how far the occupied density is from being stationary in its own field. The native convergence implementation checks both a density-change norm and a maximum-entry orbital residual, each normalized by the AO dimension. ROHF uses its effective Fock and total-density construction.

convergence_threshold defaults to 1e-7. Native dispatch uses that for the orbital-gradient condition and one hundred times it, 1e-5, for the density condition. An energy-change quantity is logged and used by optimization logic, but a small change in energy alone is not the native convergence test. Energy can appear nearly stationary while the density oscillates, particularly near shallow directions.

max_iterations=50 is the default molecular outer-iteration limit, not a scientific accuracy target. Increasing it is reasonable when residuals decrease steadily. It is not a convincing remedy for a repeating occupation cycle, an invalid guess, NaNs, or a wrong charge. Native failure is explicit; keep the error message rather than substituting a plausible energy.

The PySCF adapter’s same-named tolerance is not numerically equivalent: it sets mf.conv_tol = 0.1 * convergence_threshold, without explicitly setting mf.conv_tol_grad. Its regular run path does not itself check mf.converged before returning the energy/wavefunction pair. Consequently, a returned tuple is not a uniform convergence certificate across implementations. The adapter source, not the name of the setting, determines the contract.

Numerical consistency also depends on integral screening. In the native direct engine, a nonpositive eri_threshold triggers the rule convergence_threshold * 1e-5; a positive value is used explicitly. shell_pair_threshold defaults to 1e-12. These omit small integral contributions; neither is a guaranteed final energy-error bound. Tightening an SCF residual while retaining excessively loose integrals may request an accuracy the underlying field cannot supply.

Near-linear dependence is another issue, not an optimizer personality. The native overlap orthogonalization can remove near-redundant basis directions. Its low-level default threshold is 1e-6, and it is not a high-level qdk SCF settings key. Check the number of retained MOs and overlap warnings. Native ROHF currently requires the retained MO count to equal the AO count in its effective-Fock construction; removal of directions can make that route fail. Adding iterations cannot restore a missing linear-algebra capability.

DIIS: use the history of errors rather than guess a mixing fraction

DIIS, direct inversion in the iterative subspace, collects recent Fock matrices and their residuals. Flatten each residual matrix into a vector $e_i$; flattening simply assigns its entries a consistent order. DIIS seeks weights $c_i$ that minimize the residual of a linear combination while preserving an affine combination:

These weights are optimization coefficients, not probabilities; negative values can be valid. In a toy scalar example with residuals 2 and −1, weights one-third and two-thirds cancel the residual while summing to one. Real molecular residuals have many components, so exact cancellation is usually unavailable.

Define the residual-overlap matrix $B_{ij}=e_i^{\mathsf T}e_j$. Constrained minimization gives a small linear system for the weights and a Lagrange multiplier, a variable enforcing their sum:

The native DIIS implementation builds these inner products, rescales the overlap block by its largest entry, and solves the augmented system by column-pivoted QR. Its internal default history size is eight. When the history fills, it drops the oldest pair. If the residual-overlap block is zero it returns the current field; if the augmented system is judged too nearly dependent using its 1e-12 determinant test, it prunes old history rather than trust enormous unstable weights.

That last point matters physically. Several almost identical past errors provide little independent information. Keeping more of them does not necessarily improve the direction; it can make the extrapolation ill-conditioned. DIIS then diagonalizes its proposed field and reconstructs occupations. It is fast near many ordinary solutions but is not a line-search method requiring every proposed energy to decrease.

The low-level DIIS history and damping configuration should not be confused with available high-level settings. At this revision, do not pass invented diis_subspace_size or enable_damping keywords to the public qdk factory. Inspect the concrete schema before configuring a backend.

Level shifts stabilize the iteration without adding a new molecular force

A positive level_shift selects a modified trial field in the native DIIS implementation instead of its ordinary Pulay extrapolation. The default -1.0 leaves this disabled. The implemented expression is

This is an iteration device that separates occupied and unoccupied directions. With an orthonormal basis and a spin-resolved idempotent density, it leaves occupied levels unchanged and raises virtual levels. With the native RHF spin-summed density, $P=2\Pi_{\mathrm{occ}}$, it shifts occupied levels by $-\mu$ and virtual levels by $+\mu$; the gap is enlarged by $2\mu$. Using the actual density convention explains why the factor matters.

The physical field used for the SCF energy and residual is kept separate from the trial field. Nevertheless, a large shift can slow relaxation or steer the iteration to a different stationary solution. For final interpretation, especially orbital energies and a stability check that uses them, refine with unshifted settings on a fresh solver rather than assume any shifted trial eigenvalues are physical canonical energies.

The engine also uses incremental Fock updates after initial full builds. fock_reset_steps, default 1073741824, governs periodic rebuilding from the full density; it is not another iteration cap. More frequent rebuilding is a numerical consistency tool, not a correction for missing electron correlation.

GDM: optimize orbital rotations directly

DIIS proposes fields. Geometric direct minimization (GDM) instead asks which orbital change lowers the energy while automatically preserving orthonormality. Let $\kappa$ be a real antisymmetric matrix, meaning $\kappa^{\mathsf T}=-\kappa$. Define

The matrix exponential makes $U$ orthogonal, so if $C^{\mathsf T}SC=I$, the new coefficients satisfy the same condition. We optimize rotation parameters rather than repair broken orbital normalization after each step. For RHF the physically changing directions mix occupied and virtual orbitals; rotations wholly inside the occupied space leave the determinant unchanged up to phase. UHF has separate alpha and beta directions. ROHF needs blocks connecting doubly occupied, singly occupied, and unoccupied spaces.

Let $g$ be the derivative of the energy with respect to the independent rotation parameters. A Newton step would solve $H p=-g$, where $H$ is the orbital Hessian, or second-derivative matrix, and $p$ is a direction. Computing and storing the exact Hessian repeatedly is expensive. Quasi-Newton methods approximate its action from observed changes.

The native GDM code uses a limited-history BFGS two-loop recursion. BFGS is named for Broyden, Fletcher, Goldfarb, and Shanno. Each accepted step supplies an orbital displacement $s$ and a gradient change $y$. The approximation is updated to be consistent with the observed relation $H s\approx y$. Limited history retains these vector pairs instead of materializing a huge dense inverse Hessian.

Its starting diagonal approximation uses pseudo-canonical orbital gaps: it diagonalizes appropriate occupied and virtual blocks without changing the represented occupied subspaces, and transforms the history into that same coordinate system. For ordinary RHF/UHF, the positive denominator is built from

where $i$ and $a$ label occupied and virtual orbitals, the tildes denote the block-canonical energies, $\Delta E$ is the latest accepted energy change, and $\gamma$ is four for restricted closed-shell or two for unrestricted spin blocks. This is a preconditioner: it scales different directions so an energetically stiff rotation is not stepped as aggressively as a soft one. It is not an exact curvature measurement proving stability.

If the latest curvature product $s^{\mathsf T}y$ is too small or invalid, GDM clears its history. If its proposed direction has $p^{\mathsf T}g>0$, it points uphill and the history is discarded. Without usable history it falls back to the diagonal-preconditioned negative gradient. gdm_bfgs_history_size_limit defaults to fifty and must be at least one. More history costs memory proportional to the rotation-vector dimension and may not be more robust.

A line search makes each expensive trial count

Having a direction does not determine how far to move. A line search tests positive step lengths $\alpha$ along it, evaluating trial orbitals, densities, energies, and gradients. These trial evaluations rebuild real Coulomb/exchange and, for DFT, XC contributions; one GDM iteration can therefore cost several field evaluations.

The native Nocedal–Wright line search brackets a step then narrows the bracket. It uses the sufficient-decrease and strong-Wolfe tests. If $f(\alpha)=E(Ce^{\alpha\kappa(p)})$, they have the form

The first asks for enough decrease relative to the initial downhill slope; the second asks that the slope flatten sufficiently. A step that is tiny merely because it avoids all difficulty is not necessarily useful. These are local numerical conditions, not a guarantee of the global electronic minimum.

If a BFGS line search fails, native GDM retries using steepest descent. If that also fails, it accepts a zero rotation only when its gradient-norm condition is already sufficiently small; otherwise it aborts with an error. This is a useful distinction between “nothing needs to move” and “the algorithm cannot find a valid move.” The current GDM iteration also explicitly rejects MPI world sizes greater than one because of its trial-field participation limitation. Compiled MPI capability elsewhere does not make this particular optimization path distributed.

The native hybrid switch is precise

Public enable_gdm=False selects DIIS. Its default True selects DIIS followed, when appropriate, by GDM, not pure GDM from the first iteration. The actual switching condition is: exceed gdm_max_diis_iteration, or after more than one step have an absolute energy change below energy_thresh_diis_switch.

The defaults are fifty and 1e-3 Hartree respectively. The threshold must be positive; the internal DIIS limit is clamped to at least two. Once the switch occurs, the algorithm stays in GDM. Some documentation/comments describe the opposite inequality; the current condition uses <. We can understand the intended division of labor without that misleading wording: DIIS approaches a candidate solution, then direct orbital minimization refines it. Neither step proves that symmetry-forbidden downhill directions have been explored.

Stability is curvature, not another convergence tolerance

Near stationary orbitals, express the energy in independent rotation coordinates $x$:

$H_{\mathrm{orb}}$ is the electronic orbital Hessian. It is not the nuclear Hessian of Lesson 7. SCF stationarity makes $g$ small. If $H_{\mathrm{orb}}v=\lambda v$ with a negative $\lambda$, then a small motion in direction $v$ lowers the energy to second order even though the first derivative was zero.

For an illustrative toy surface, $E(u,v)=E_0+0.2u^2-0.05v^2$, the point $u=v=0$ has zero gradient but Hessian eigenvalues 0.4 and −0.1. It is a saddle. Here $u,v$ are orbital-rotation coordinates, not changes in the N–N bond. No numerical conclusion about our nitrogen calculation follows from these invented coefficients.

An internal instability lowers energy inside the same orbital restriction, such as RHF to a different RHF determinant. An external instability relaxes a restriction, here typically common alpha/beta orbitals becoming different: RHF to UHF or the corresponding KS channel. Native external checking is not an exhaustive search over all complex, noncollinear, or spinor states. The PySCF external helper explicitly excludes real-to-complex instability.

Even an internally and externally stable result is only locally stable in the checked channels and within the selected numerical thresholds. It can still be above another minimum far away. A lower unrestricted solution may be spin-contaminated; it must not simply be relabeled a spin-pure singlet because it has seven electrons of each spin.

What the stability checkers actually compute

The interface is create("stability_checker", implementation).run(wavefunction), returning (is_stable, StabilityResult). The native implementation is "qdk"; the optional plugin is "pyscf".

The native checker uses orbital energies, coefficients, a basis, and electron counts to build an occupied–virtual response problem. It reconstructs Aufbau densities, so use this for an appropriate SCF reference, not as a stability theory for arbitrary correlated coefficients that happen to fit in Wavefunction. It does not replace the converged reference requirement with an independent stationarity certification.

It applies the Hessian to trial vectors without building the entire Hessian matrix. A trial orbital variation induces a trial density; Coulomb/exchange and, for DFT, the XC response kernel give its field response. Adding orbital-energy-gap terms completes the action. Davidson diagonalization iteratively finds low eigenvalues in a small expanding subspace, preconditioned by diagonal information. This native eigensolver is supplied by MACIS; the integral and XC-response libraries provide the expensive operator actions. These are concrete library boundaries, not separate chemical approximations.

Native defaults are internal=True, external=False, max_subspace=80, davidson_tolerance=1e-8, stability_tolerance=-1e-4, and method="hf". It finds one lowest eigenpair per requested channel; max_subspace limits search-space size and is clamped to the rotation-problem dimension. It is not a setting for how many physical unstable states to report. Davidson failure is an eigensolver failure, not evidence of stability; a larger useful subspace or better-conditioned reference may be needed.

For our unreduced RHF model, seven occupied and 21 virtual orbitals give $7\times21=147$ independent occupied–virtual rotation coordinates. That source-derived dimension is not a computed eigenvalue. UHF concatenates alpha and beta rotation blocks. Returned vectors use the public occupied-index-fast convention; use the provided rotation utility instead of guessing how to reshape a raw vector.

The classification is lowest_eigenvalue > stability_tolerance. With the default negative threshold, a weakly negative value such as −0.00005 passes this numerical criterion, while −0.0002 does not. These are illustrative values. The Davidson solve tolerance and the physical decision threshold serve different purposes and should be assessed together near zero curvature.

Native checking rejects ROHF/ROKS references even though native SCF can produce ROHF. It also rejects external checking on unrestricted references. Requesting neither channel is an error. Its explicit zero-electron case returns stable status with empty eigenpair arrays. For other results, inspect has_internal_result() and has_external_result() before requesting eigenvalues; an unchecked channel has a default status, not a completed calculation. StabilityResult provides the two statuses, eigenvalue arrays, eigenvector columns, and smallest-eigenpair getters. It is a result report, not proof that every possible direction was checked.

The PySCF checker explicitly requires a single-determinant state-vector container. It supports internal RHF, ROHF and UHF analysis and RHF external analysis, including corresponding supported KS response paths. External ROHF/UHF requests raise errors. Its defaults include external=True, nroots=3, with_symmetry=False, davidson_tolerance=1e-8, stability_tolerance=-1e-4, method="hf", xc_grid=3, and pyscf_verbose=4. with_symmetry controls point-group restrictions where supplied by PySCF, not whether alpha and beta are identical.

This adapter rescales internal RHF eigenvalues by one-quarter and UHF eigenvalues by one-half to align reporting conventions; external values are left unscaled. Comparing raw eigenvalues from different implementations without matching conventions is therefore unsafe. In contrast to native checking, disabling both PySCF channels returns empty data with stable defaults. That documented-by-test behavior means “nothing requested,” not “minimum established.” Always match the checker’s method and relevant grid settings to the SCF model; it does not infer the functional reliably from a generic wavefunction object.

Perturb the orbitals, then solve again

A negative eigenvector is a direction, not a new converged state. The rotation utility unpacks the vector into an antisymmetric matrix, exponentiates it, and returns new Orbitals. It invalidates orbital energies because rotated orbitals are no longer guaranteed canonical.

For an RHF external perturbation, restricted_external=True constructs unrestricted orbitals by rotating alpha while retaining the original beta coefficients. Internal restricted perturbations rotate their common orbital set. For unrestricted input the vector contains separate alpha and beta blocks. The input assumes Aufbau occupations.

Here is a complete single-retry illustration using reference, orbitals, and structure from Lesson 5. The factor 0.1 is an explicitly illustrative perturbation amplitude, not a native setting or a universally good step:

from qdk_chemistry.algorithms import create
from qdk_chemistry.utils import rotate_orbitals

checker = create("stability_checker", "qdk", method="hf", external=True)
stable, report = checker.run(reference)
if not stable:
    external = report.has_external_result() and not report.is_external_stable()
    if external:
        eigenvalue, direction = report.get_smallest_external_eigenvalue_and_vector()
    else:
        eigenvalue, direction = report.get_smallest_internal_eigenvalue_and_vector()
    n_alpha, n_beta = reference.get_total_num_electrons()
    guess = rotate_orbitals(orbitals, 0.1 * direction, n_alpha, n_beta, external)
    retry = create(
        "scf_solver", "qdk", method="hf",
        scf_type="unrestricted" if external else "restricted",
    )
    retry_energy, retry_reference = retry.run(structure, 0, 1, guess)
    recheck = create("stability_checker", "qdk", method="hf", external=not external)
    retry_stable, retry_report = recheck.run(retry_reference)

The recheck matters: neither the perturbation nor a lower returned energy proves final stability. Changing the perturbation scale affects which basin the reoptimization reaches. A very large rotation need not lower the energy predicted by a local curvature expansion. New algorithm instances are used when their settings change; run locks ordinary settings before execution, including failing execution. As Lesson 3 explained, that is not a guarantee of deep freezing every separately reachable nested object.

Automated stabilized SCF has a specific policy

create("scf_solver", "qdk_stabilized") performs the repeated workflow. Its implementation holds nested AlgorithmRef settings for scf_solver and stability_checker, both native by default. The defaults are five stability cycles, both check flags true, fail_on_unstable=True, and external_instability_action="unrestricted".

The native wrapper runs SCF, checks applicable channels, prefers an external unstable direction when both are present, rotates, and reruns SCF. Unlike our manual illustration, it passes the returned direction without an additional exposed amplitude multiplier. When it follows an external direction, its default policy forces subsequent SCF to unrestricted. Once the reference is unrestricted, the wrapper stops requesting the unsupported RHF external channel and checks internal stability.

max_stability_iterations has a declared nonnegative bound. With a positive value $m$, at most $m$ rotation/rerun cycles follow the initial SCF; a final check evaluates the last result. With zero, the wrapper is plain SCF passthrough and does no stability check. fail_on_unstable=False permits return of the last energy/wavefunction even if the requested local stability was not achieved. The return type remains the ordinary SCF pair, not a status-rich stability report; run a final explicit checker if that report is needed.

The alternative action "rotate_only" suppresses the wrapper’s explicit switch to unrestricted. It does not make the generated external guess restricted again. A native auto/restricted singlet child can reject that unrestricted guess, so this is not a safe synonym for “remove external instability while remaining RHF.” Also, the default native checker cannot analyze a ROHF result; select a compatible checker or a manual supported workflow rather than assume the composite expands its components’ capabilities.

A concrete configuration example, still unexecuted, is:

from qdk_chemistry.data import AlgorithmRef

robust = create(
    "scf_solver", "qdk_stabilized",
    method="hf", max_iterations=100, max_stability_iterations=3,
    fail_on_unstable=True,
    scf_solver=AlgorithmRef("scf_solver", "qdk", enable_gdm=True),
    stability_checker=AlgorithmRef(
        "stability_checker", "qdk", max_subspace=80,
        davidson_tolerance=1e-8, stability_tolerance=-1e-4,
    ),
)
robust_energy, robust_reference = robust.run(structure, 0, 1, "cc-pvdz")

Common parent SCF settings are copied to the child after nested construction, so the parent’s max_iterations wins over a competing nested override. The parent also supplies check flags and method to a compatible checker. Disabling both checks with positive cycles is not meaningful and will cause the native checker to reject the request.

The optional "pyscf_stabilized" implementation is different composition, not an alias for this native wrapper. It calls the PySCF mean-field object’s own stability routine, obtains candidate orbital rotations/status, preferentially follows an available external instability, converts to an unrestricted mean-field object when necessary, and reruns kernel with the new density. It does not call the separately registered QDK "pyscf" stability checker or inherit all that checker’s nroots/threshold controls. Its defaults are five cycles, both checks enabled, and failure on unresolved instability; zero cycles is passthrough. Its cycle setting is an integer without the native wrapper’s declared nonnegative bound, so supply a sensible nonnegative value explicitly. The regular adapter’s missing explicit mf.converged check is not repaired merely by naming the wrapper “stabilized.”

A practical diagnostic record for stretched nitrogen

For our 1.85-Å case, first retain the exact geometry, electron count, basis definition, SCF type, method, backend, and configured thresholds. Check the iteration log and retained basis rank before interpreting the last decimal places. If residuals decrease steadily, a larger iteration cap may be justified. If they oscillate, compare the default hybrid with a deliberate alternative or a compatible orbital guess; do not change basis, functional, and optimizer simultaneously and then call the difference a convergence correction.

Next verify numerical consistency by tightening relevant thresholds and comparing converged energies and density-sensitive quantities, rather than assuming a requested tolerance is an error bar. Examine internal stability of a properly converged reference and external stability if the chosen reference and scientific question permit it. Keep the original restricted result when following an unrestricted branch: this is a change in allowed variational space and potentially spin character, not just deletion of a failed file.

The repository’s stability tests include RHF nitrogen at 1.2 and 1.6 Å in def2-SVP with different internal/external status patterns. The stabilized-SCF test checks external-priority reoptimization at its own 1.6-Å def2-SVP model. Those are recorded source baselines establishing that these code paths matter; they are not the missing results of our cc-pVDZ 1.85-Å calculation. The standalone documentation workflow uses yet another geometry, 1.4 Å, and is a procedural example rather than our reference.

Finally, report exactly what was established: numerical convergence, the channels and curvature tolerance checked, whether reoptimization changed spin restrictions, and any remaining unresolved status. A stable HF reference remains an approximation. We now have a disciplined starting point for interpreting density and for asking a different question: how the energy changes when the nuclei, rather than the orbitals alone, move.