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.

Ask a Question

Have a question about quantum chemistry, or anything discussed on this site? Send it our way and we will reply by email.

Prefer to fill it out directly? Open the form.