4.5.9 Internal Stability Analysis and Automated Correction for Energy Minima

At convergence, the SCF energy will be at a stationary point with respect to changes in the MO coefficients. However, this stationary point is not guaranteed to be an energy minimum, and in cases where it is not, the wave function is said to be unstable. Even if the wave function is at a minimum, this minimum may be an artifact of the constraints placed on the form of the wave function. For example, an unrestricted calculation will usually give a lower energy than the corresponding restricted calculation, and this can give rise to an RHF $\to$ UHF instability.

Based on our experience, even for very simple data set such as the G2 atomization energies,²⁰³ using the default algorithm (DIIS) produces unstable solutions for several species (even for single atoms with some density functionals). In such cases, failure to check the internal stability of SCF solutions can result in flawed benchmark results. Although in general the use of gradient-based algorithms such as GDM is more likely to locate the true minimum, it still cannot entirely eliminate the possibility of finding an unstable solution.

To understand what instabilities can occur, it is useful to consider the most general form possible for the spin orbitals:

$${\chi }_i(𝐫,\zeta )={\psi }_i^{\alpha }(𝐫)\alpha (\zeta )+{\psi }_i^{\beta }(𝐫)\beta (\zeta ).$$

(4.39)

Here, ${\psi }_i^{\alpha }$ and ${\psi }_i^{\beta }$ are complex-valued functions of the Cartesian coordinates $𝐫$, and $\alpha$ and $\beta$ are spin eigenfunctions of the spin-variable $\zeta$. The first constraint that is almost universally applied is to assume the spin orbitals depend only on one or other of the spin-functions $\alpha$ or $\beta$. Thus, the spin-functions take the form

$${\chi }_i(𝐫,\zeta )={\psi }_i^{\alpha }(𝐫)\alpha (\zeta )\mathit{\hspace{1em}}\mathrm{or}\mathit{\hspace{1em}}{\chi }_i(𝐫,\zeta )={\psi }_i^{\beta }(𝐫)\beta (\zeta ).$$

(4.40)

In addition, most SCF calculations use real functions, and this places an additional constraint on the form of the wave function. If there exists a complex solution to the SCF equations that has a lower energy, the wave function exhibits a real $\to$ complex instability. The final constraint that is commonly placed on the spin-functions is that ${\psi }_i^{\alpha }={\psi }_i^{\beta }$, i.e., that the spatial parts of the spin-up and spin-down orbitals are the same. This gives the familiar restricted formalism and can lead to an RHF $\to$ UHF instability as mentioned above. Further details about the possible instabilities can be found in Ref. 829.

Wave function instabilities can arise for several reasons, but frequently occur if

• •

  There exists a singlet diradical at a lower energy then the closed-shell singlet state.

• •

  There exists a triplet state at a lower energy than the lowest singlet state.

• •

  There are multiple solutions to the SCF equations, and the calculation has not found the lowest energy solution.

Q-Chem’s previous stability analysis package suffered from the following limitations:

• •

  It is only available for restricted (close-shell) and unrestricted SCF calculations.

• •

  It requires the analytical orbital Hessian of the wave function energy.

• •

  The calculation terminates after the corrected MOs are generated, and a second job is needed to read in these orbitals and run another SCF calculation.

The implementation of internal stability analysis in GEN_SCFMAN overcomes almost all these shortcomings. Its availability has been extended to all the implemented orbital types. As in the old code, when the analytical Hessian of the given orbital type and theory (e.g. RO/B3LYP) is available, it computes matrix-vector products analytically for the Davidson algorithm.²¹³ If the analytical Hessian is not available, users can still run stability analysis by using the finite-difference matrix-vector product technique developed by Sharada et al.,⁸⁴¹ which requires the gradient (related to the Fock matrix) only:

$${\mathrm{𝐇𝐛}}_1=\frac{\nabla E({𝐗}_0+\xi {𝐛}_1)-\nabla E({𝐗}_0-\xi {𝐛}_1)}{2\xi }$$

(4.41)

where $𝐇$ is the Hessian matrix, ${𝐛}_1$ is a trial vector, ${𝐗}_0$ stands for the current stationary point, and $\xi$ is the finite step size. With this method, internal stability analysis is available for all the implemented orbital types in GEN_SCFMAN. It should be noted that since the second derivative of NLC functionals such as VV10 is not available in Q-Chem, this finite-difference method will be used by default for the evaluation of Hessian-vector products.

GEN_SCFMAN allows multiple SCF calculations and stability analyses to be performed in a single job so that it can make use of the corrected MOs and locate the true minimum automatically. The MOs are displaced along the direction of the lowest-energy eigenvector (with line search) if an SCF solution is found to be unstable. A new SCF calculation that reads in these corrected MOs as initial guess will be launched automatically if INTERNAL_STABILITY_ITER > 0. Such macro-loops will keep going until a stable solution is reached.

Note:  The stability analysis package can be used to analyze both HF and DFT wave functions.

$molecule
   0 3
   b
   b 1 R

   R = 1.587553
$end

$rem
   JOBTYPE              sp
   METHOD               b97m-v
   BASIS                6-31g
   UNRESTRICTED         true
   THRESH               14
   SYMMETRY             false
   SYM_IGNORE           true
   SCF_FINAL_PRINT      1
   SCF_ALGORITHM        gdm
   SCF_CONVERGENCE      8
   INTERNAL_STABILITY   true  !turn on internal stability analysis
   FD_MAT_VEC_PROD      false !use finite-diff for the vv10 part only
$end