Stability of Wavefunctions

Wavefunctions generated by SCF calculations can be unstable in various ways:

1) The lowest energy wavefunction is a singlet biradical instead of a closed shell singlet. A proper description of singlet biradicals at the Hartree-Fock level requires an UHF wavefunction. This is a typical case for an RHF/UHF instability

2) A triplet state is energetically more favorable than the lowest lying singlet state. This will also lead to an RHF/UHF instability

3) There is more than one solution to the SCF equations for the system and the calculation converges to a less favorable one. This will lead to either a RHF/RHF or UHF/UHF instability

The various cases of wavefunction instability will be demonstrated using two small model systems, O₂ and O₃.

1) molecular oxygen O₂
Let us first calculate the Hartree-Fock energy for singlet dioxygen at the RHF/STO-5G level of theory (using the geometry optimized at this level) with the following input file:

#RHF/STO-5G scf=(tight,qc)

RHF/STO-3G singlet O2

0 1
O1
O2  1  r1

r1=1.22220

O=O

This yields the following result:

 SCF Done:  E(RHF) =  -148.886061396     a.u. after    4 cycles
            Convg  =    0.1803D-07                    16 Fock formations.
              S**2 =  0.0000                  -V/T =  2.0033

 **********************************************************************

            Population analysis using the SCF density.

 **********************************************************************

 Orbital symmetries:
       Occupied  (SGG) (SGU) (SGG) (SGU) (PIU) (SGG) (PIU) (PIG)
       Virtual   (PIG) (SGU)
 Unable to determine electronic state:  partially filled degenerate orbitals.

Subsequent analysis of this wavefunction with the keyword combination stable guess=read reads in the converged wavefunction and, by analysis of a number of excitations starting from the HF-solution, yields the following additional information:

 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Excited state symmetry could not be determined.
 Eigenvector   1:   Triplet-?Sym  Eigenvalue=-0.2009390
       7 ->  9         0.70536

 Excited state symmetry could not be determined.
 Eigenvector   2:   Triplet-?Sym  Eigenvalue=-0.1078684
       8 ->  9         0.70711

 Excited state symmetry could not be determined.
 Eigenvector   3:   Singlet-?Sym  Eigenvalue= 0.0000000
       8 ->  9         0.70711

 Excited state symmetry could not be determined.
 Eigenvector   4:   Triplet-?Sym  Eigenvalue= 0.0667651
       6 ->  9         0.65552
       7 -> 10         0.26170

 Excited state symmetry could not be determined.
 Eigenvector   5:   Triplet-?Sym  Eigenvalue= 0.2052332
       5 ->  9         0.70711

 Excited state symmetry could not be determined.
 Eigenvector   6:   Singlet-?Sym  Eigenvalue= 0.2170234
       6 ->  9         0.66826
       7 -> 10         0.22006
 The wavefunction has an RHF -> UHF instability.

This analysis points to a triplet state wavefunction lower in energy than the current singlet state. The actual triplet wavefunction is, however, not calculated explicitly. In order to find the optimized triplet wavefunction, a second calculation must be performed. Using the same geometry as before (RHF/STO-5G), the calculation is performed in one run together with the stability analysis of the triplet wavefunction:

#UHF/STO-5G stable scf=(tight,qc)

UHF/STO-3G triplet O2

0 3
O1
O2  1  r1

r1=1.22220

The total energy of -148.968737 Hartree obtained in this UHF/STO-5G calculation is lower by 217 kJ/mol than the one obtained for the singlet state at the RHF/STO-3G level! Despite this large improvement, the stability analysis reveals one more problem with the wavefunction:

 The wavefunction has an internal instability.

From the additional data given by the stability analysis it appears that the triplet state optimized with the default guess does NOT converge to the triplet state of correct symmetry (and lowest energy). This can be solved either by an appropriate manipulation of the initial guess as discussed before or with the keyword guess=mix originally designed to provide an asymmetric initial guess for calculations on singlet biradicals. In either case, a new triplet state of different symmetry is obtained at somewhat lower total energy of -148.9720434, 225.7 kJ/mol below the singlet state. Stability analysis of this wavefunction now confirms that:

 The wavefunction is stable under the perturbations considered.

The same result can be obtained by running the stability calculation with stable=opt. Under these latter conditions the constraints imposed upon the wavefunction are reduced incrementally until a stable wavefunction is obtained. The key feature of the most stable wavefunction obtained for triplet oxygen by either stable=opt or guess=mix is the reduced symmetry of the two highest lying orbitals.

2) ozone O₃
Calculation of the singlet state RHF/6-31G(d) energy of ozone at its experimental geometry with the following input

#RHF/6-31G(d) stable scf=(tight,qc)

RHF/STO-3G singlet O3

0 1
O1
O2  1  1.2227
O3  1  1.2227  2  114.0451

and subsequent stability analysis yields the following result:

 SCF Done:  E(RHF) =  -224.258537798     a.u. after    8 cycles
 The wavefunction has an RHF -> UHF instability.

As we have seen in the first example of O₂, this can immediately be solved by performing an UHF/6-31G(d) calculation on the corresponding triplet state of ozone. In contrast to the first example, however, the total energy of the triplet state obtained in this way is higher than that of the singlet state obtained initially. Stability analysis also reveals an internal instability as encountered before due to convergence to a state of wrong symmetry:

 SCF Done:  E(UHF) =  -224.226611298     a.u. after   10 cycles
 The wavefunction has an internal instability.

In this situation one must consider the possibility of a singlet diradical which requires the use of an UHF wavefunction even for a singlet state. Generating a broken symmetry guess for the singlet wavefunction with guess=mix currently appears to work well only with the INDO guess:

#UHF/6-31G(d) guess=(INDO,mix) stable scf=(tight,qc)

UHF/STO-3G singlet O3

0 1
O1
O2  1  1.2227
O3  1  1.2227  2  114.0451

As a result a biradical singlet state is obtained that is 180.3 kJ/mol more favorable than the singlet state described by the RHF wavefunction before. In this last case the stability analysis detects no further problems with the wavefunction:

SCF Done:  E(UHF) =  -224.327207261     a.u. after   10 cycles
The wavefunction is stable under the perturbations considered.

Please observe that a broken symmetry UHF singlet wavefunction is only obtained using the guess=mix keyword. If this is not used, the initial guess is chosen such that the SCF calculation converges to the RHF wavefunction even with the UHF Ansatz. An alternative strategy to obtain the energetically most favorable singlet wavefunction for ozone involves the use of stable=opt.