Skip to content

Davidson Iterative Diagonalization

For a molecule with \(n_{occ}\) occupied and \(n_{virt}\) virtual orbitals, the stability matrix has dimension \(n_{occ} \times n_{virt}\), which grows quickly into the thousands even for modest molecules. Forming and diagonalizing it directly is computationally impossible for routine use — Gaussian instead uses an iterative Davidson-type algorithm (the CISAX routine seen in the output) to extract only the lowest few eigenvalues, which are all that stability actually requires.

The procedure, in outline:

  1. Initial trial vectors. Start from a small set of guess vectors \(b_1, b_2, \dots\) spanning a subspace \(V\) — typically built from the orbital-energy-ordered occupied → virtual pairs.
  2. Matrix–vector products. Instead of forming the full Hessian, compute only the action of the Hessian on each trial vector, \( \sigma_k = H b_k \), using the two-electron integrals directly.
  3. Reduced (projected) matrix. Build and diagonalize the small matrix \( \tilde{H}_{kl} = b_k^{\dagger} H b_l \) within the current subspace — cheap, since its dimension equals the number of trial vectors, not the full problem size.
  4. Residuals. For each approximate eigenpair \((\lambda, v)\) from the reduced matrix, compute the residual \( r = Hv - \lambda v \). If \( \|r\| \) is below the convergence threshold (0.001 in the worked example), that root is converged.
  5. Subspace expansion. For unconverged roots, form a new trial vector from the residual (preconditioned by the orbital energy differences) and add it to \(V\).
  6. Repeat steps 2–5 until all requested roots converge — visible in the log as the Dimension growing from 24 to 61 across five iterations for the formic acid dimer example.

This iterative strategy is what lets Gaussian test stability for systems where the full Hessian would be far too large to handle directly.