X

Search Results

Searching....

4.5 Converging SCF Calculations

4.5.16 Internal Stability Analysis and Automated Correction for Energy Minima

(April 13, 2024)

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 UHF instability.

Based on our experience, even for very simple data set such as the G2 atomization energies, 262 Curtiss L. A. et al.
J. Chem. Phys.
(1997), 106, pp. 1063.
Link
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:

χi(𝐫,ζ)=ψiα(𝐫)α(ζ)+ψiβ(𝐫)β(ζ). (4.45)

Here, ψiα and ψiβ are complex-valued functions of the Cartesian coordinates 𝐫, and α and β are spin eigenfunctions of the spin-variable ζ. The first constraint that is almost universally applied is to assume the spin orbitals depend only on one or other of the spin-functions α or β. Thus, the spin-functions take the form

χi(𝐫,ζ)=ψiα(𝐫)α(ζ)orχi(𝐫,ζ)=ψiβ(𝐫)β(ζ). (4.46)

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 complex instability. The final constraint that is commonly placed on the spin-functions is that ψiα=ψiβ, 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 UHF instability as mentioned above. Further details about the possible instabilities can be found in Ref.  1105 Seeger R., Pople J. A.
J. Chem. Phys.
(1977), 66, pp. 3045.
Link
.

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. 275 Davidson E. R.
J. Comput. Phys.
(1975), 17, pp. 87.
Link
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., 1122 Sharada S. M. et al.
Mol. Phys.
(2015), 113, pp. 1802.
Link
which requires the gradient (related to the Fock matrix) only:

𝐇𝐛1=E(𝐗0+ξ𝐛1)-E(𝐗0-ξ𝐛1)2ξ (4.47)

where 𝐇 is the Hessian matrix, 𝐛1 is a trial vector, 𝐗0 stands for the current stationary point, and ξ 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 (except for 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.

4.5.16.1 Job Control

INTERNAL_STABILITY

INTERNAL_STABILITY
       Perform internal stability analysis in GEN_SCFMAN.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not perform internal stability analysis after convergence. TRUE Perform internal stability analysis and generate the corrected MOs.
RECOMMENDATION:
       Turn it on when the SCF solution is prone to unstable solutions, especially for open-shell species.

FD_MAT_VEC_PROD

FD_MAT_VEC_PROD
       Compute Hessian-vector product using the finite difference technique.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE (TRUE when the employed functional contains non-local correlation (except VV10))
OPTIONS:
       FALSE Compute Hessian-vector product analytically. TRUE Use finite difference to compute Hessian-vector product.
RECOMMENDATION:
       Set it to TRUE when analytical Hessian is not available. Note:  For simple R and USCF calculations, it can always be set to FALSE, which indicates that only the NLC part will be computed with finite difference (if its analytic orbital hessian is unavailable).

INTERNAL_STABILITY_ITER

INTERNAL_STABILITY_ITER
       Maximum number of new SCF calculations permitted after the first stability analysis is performed.
TYPE:
       INTEGER
DEFAULT:
       0 (automatically set to 1 if INTERNAL_STABILITY = TRUE)
OPTIONS:
       n n new SCF calculations permitted.
RECOMMENDATION:
       Give a larger number if 1 is not enough (still unstable).

INTERNAL_STABILITY_DAVIDSON_ITER

INTERNAL_STABILITY_DAVIDSON_ITER
       Maximum number of Davidson iterations allowed in one stability analysis.
TYPE:
       INTEGER
DEFAULT:
       50
OPTIONS:
       n Perform up to n Davidson iterations.
RECOMMENDATION:
       Use the default.

INTERNAL_STABILITY_CONV

INTERNAL_STABILITY_CONV
       Convergence criterion for the Davidson solver (for the lowest eigenvalues).
TYPE:
       INTEGER
DEFAULT:
       4 (3 when FD_MAT_VEC_PROD = TRUE)
OPTIONS:
       n Terminate Davidson iterations when the norm of the residual vector is below 10-n.
RECOMMENDATION:
       Use the default.

INTERNAL_STABILITY_ROOTS

INTERNAL_STABILITY_ROOTS
       Number of lowest Hessian eigenvalues to solve for.
TYPE:
       INTEGER
DEFAULT:
       2
OPTIONS:
       n Solve for n lowest eigenvalues.
RECOMMENDATION:
       Use the default.

Example 4.18  Unrestricted SCF calculation of triplet B2 using B97M-V/6-31g with the GDM algorithm. A displacement is performed when the first solution is characterized as a saddle point, and the second SCF gives a stable solution.

$molecule
   0 3
   b
   b 1 R

   R = 1.587553
$end

$rem
   JOBTYPE              sp
   METHOD               b97m-v
   BASIS                6-31g
   UNRESTRICTED         true
   THRESH               14
integral_symmetry false
point_group_symmetry False
   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