Skip to content

How Does a Stability Analysis Work?

An SCF calculation stops as soon as the energy is stationary with respect to mixing occupied and virtual orbitals — the Hartree–Fock (or Kohn–Sham) analogue of Brillouin's theorem. That condition, however, only guarantees a stationary point in orbital-rotation space. It does not guarantee that the point is a minimum.

A converged determinant can fall into one of several categories:

  • Genuinely stable — a true local minimum; no lower-energy rearrangement of the orbitals exists.
  • RHF→RHF unstable — a lower-energy solution exists within the same spin-restricted, real-orbital space, often caused by an artificial symmetry constraint.
  • RHF→UHF unstable — a lower-energy spin-unrestricted solution exists; common for biradicaloid systems, stretched bonds, and some transition-metal complexes.
  • Complex unstable — a lower-energy solution exists only when complex-valued orbital coefficients are allowed; rare, and mostly relevant to certain symmetric or periodic systems.

The Orbital Hessian

To tell these cases apart, Gaussian examines the second derivative of the energy with respect to orbital-rotation parameters — the electronic (orbital) Hessian, commonly called the stability matrix. This object is built from exactly the same two blocks, A and B, that appear in RPA/TDHF and CIS excited-state theory, which is why stability analysis and excited-state calculations share so much underlying machinery.

  • If the relevant combination of A and B is positive definite, the tested channel is stable.
  • If it has a negative eigenvalue, that eigenvector points directly toward a lower-energy determinant.

Why Iterative Diagonalization?

The full stability matrix has a dimension of (occupied × virtual orbitals), which becomes far too large to diagonalize directly for any realistic molecule. Gaussian instead uses a Davidson-type iterative solver (visible in the log as the CISAX routine): a small trial subspace is built, its projected matrix is diagonalized, and new trial vectors are added from the residuals until only the lowest few eigenvalues — the ones that matter for stability — have converged.

Converged SCF Wavefunction
        │
        ▼
Form Trial Orbital-Rotation Vectors
        │
        ▼
Davidson Iterations:
  build reduced matrix → diagonalize →
  check residuals → expand subspace
        │
        ▼
Lowest Eigenvalues Converged
        │
        ▼
 Classify by Spin (⟨S²⟩) and Symmetry
        │
        ▼
 All Eigenvalues ≥ 0 ?
      │            │
     Yes            No
      │            │
   Stable      Unstable
      │            │
    Done      Rotate orbitals along
              softest eigenvector,
              restart SCF, repeat

What "Stable=Opt" Adds

Used alone, the stable keyword defaults to Stable=Opt. If any instability is found, Gaussian does not just report it — it automatically displaces the orbitals along the softest (most negative) eigenvector, forms a new SCF guess from the rotated orbitals, and reconverges. This whole cycle repeats until a genuinely stable wavefunction is found or an internal iteration limit is reached, so that the wavefunction handed off to any subsequent calculation is trustworthy.