Is an SCF solution a minimum? Orbital-rotation stability

An SCF run stops where the density commutes with its own Fock matrix. That condition says the energy is stationary: no small rotation of an occupied orbital into a virtual one changes E to first order. It does not say E goes up in every direction. Stretched water in 6-31G has an RHF solution, the one PySCF reaches from its default guess, that is a saddle point: rotating the orbitals the right way lowers E by 0.028 Eh, down to a second RHF solution. Triplet O2 in STO-3G has the same kind of UHF saddle, 1.3 mEh above the minimum. The iteration cannot tell these states apart, since both are fixed points of it, so this module asks the second-order question directly.

The rotation. Take orthonormal orbitals C (C^T S C = I) and an antisymmetric matrix K whose only nonzero entries couple an occupied orbital i to a virtual one a: K_ai = kappa_ai, K_ia = -kappa_ai. The rotated orbitals C exp(K) stay orthonormal for every kappa, so E(kappa) is the energy of a legal determinant everywhere, and

E(kappa) = E(0) + g . kappa + (1/2) kappa^T Hess kappa + O(kappa^3).

A state is stable when every eigenvalue of Hess is above -STAB_TOL (10⁻⁶). The lowest eigenvector of an unstable state is the direction in which E falls fastest.

The Hessian, RHF (real to real). The density is D = 2 C_occ C_occ^T. Expanding exp(K) to second order and E(D) to second order in the density change gives, with MO indices i, j occupied and a, b virtual and F the Fock matrix in the MO basis,

g_ai        = 4 F_ai
Hess_ai,bj  = 4 [ F_ab d_ij - F_ij d_ab + 4 (ai|bj) - (ab|ij) - (aj|bi) ].

The F_ab d_ij - F_ij d_ab part comes from the second-order change of the density along exp(K) (it moves weight from occupied to virtual orbitals); the integrals come from the Coulomb and exchange response to the first-order change. Because the F terms carry the second-order part of exp(K) exactly, this is the second derivative along the path C exp(t K) even where the gradient is not zero, and a finite difference along that path can check it anywhere, not only at a converged state.

The Hessian, UHF. Each spin has its own rotation kappa^s and density D_s = C_s,occ C_s,occ^T, so every factor of 2 from double occupation goes:

g^s_ai            = 2 F^s_ai
Hess^ss_ai,bj     = 2 (F^s_ab d_ij - F^s_ij d_ab) + 4 (ai|bj) - 2 (ab|ij) - 2 (aj|bi)
Hess^ab_ai,bj     = 4 (a_a i_a | b_b j_b)

The alpha-beta block has Coulomb only, because exchange acts within one spin. Setting kappa^a = kappa^b on a closed shell gives back the RHF energy, so the RHF Hessian is the sum of the four UHF blocks; a test checks this.

Two ways to build it. stab_hess_mo transforms the integrals to the MO basis and fills the formulas above. stab_hess_jk never forms MO integrals: it applies the Hessian to each unit vector through a J/K build on the symmetric AO matrix X = C (kappa + kappa^T) C^T,

(Hess kappa)_ai = 4 (F_vv kappa - kappa F_oo)_ai + 8 [C^T (J(X) - K(X)/2) C]_ai          RHF
(Hess kappa)^s_ai = 2 (F_vv kappa - kappa F_oo)^s_ai + 2 [C_s^T (J(X_a + X_b) - K(X_s)) C_s]_ai   UHF

which is how large codes apply it (PySCF's gen_response). The two share nothing but the Fock matrix and must agree entry for entry; a third check, the finite-difference curvature of E(C exp(t K)) along the lowest eigenvector, involves neither.

Follow-down. From an unstable state, stab_follow rotates the orbitals along the lowest eigenvector, picks the step on a grid of STAB_NSTEP points per sign (spacing STAB_DT radians) that gives the lowest energy, and restarts the SCF from the rotated density. It repeats until the lowest eigenvalue is above -STAB_TOL or STAB_MAXFOLLOW restarts have run. The first state gets the label STAB_DEFAULT (0) when it is a saddle and STAB_STABLE (1) when it is already a minimum; the end state gets its own label.

Result of stab_follow(p, D0) and stab_run(p, r), a Stability record:

lowest       lowest Hessian eigenvalue of the input state
label        its label: STAB_DEFAULT or STAB_STABLE
energy       E after follow-down (the input E when stable)
label_end    label of the end state
lowest_end   its lowest eigenvalue
steps        SCF restarts taken
energy0      E of the input state
converged    1 when every restart converged
density      the end densities, [nd, n, n]

Orbital pack. The orbitals travel as orb[nd, n + 1, m]: rows 0..n-1 of slab d are the MO coefficients of spin d (or the closed-shell orbitals when nd = 1), and row n holds 1 for an occupied orbital and 0 for a virtual one, so a Hessian member needs no other input to know the occupations. Occupations are aufbau from F(D), which is what the SCF driver uses.

Both Hessian members have the signature (orb: Arr, H: Arr, Gp: Arr) -> [np, np] and call jk_packed directly, so like the guess gs_avg they assume in-core packed integrals. ScfOptions.stab_hessian picks the member the follow-down uses. Everything here takes the run's ScfSetup (lib/hf/scf_loop) for the matrices, the electron counts and the options.

Used with import hf.stab.