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 HartreeFock 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_{2} and O_{3}.
1) molecular oxygen O_{2}
Let us first calculate the HartreeFock energy for singlet dioxygen at the RHF/STO5G level of theory (using the geometry optimized at this level) with the following input file:
#RHF/STO5G scf=(tight,qc)
RHF/STO3G 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.1803D07 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 HFsolution, 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/STO5G), the calculation is performed in one run together with the stability analysis of the triplet wavefunction:
#UHF/STO5G stable scf=(tight,qc)
UHF/STO3G triplet O2
0 3
O1
O2 1 r1
r1=1.22220
The total energy of 148.968737 Hartree obtained in this UHF/STO5G calculation is lower by 217 kJ/mol than the one obtained for the singlet state at the RHF/STO3G 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_{3}
Calculation of the singlet state RHF/631G(d) energy of ozone at its experimental geometry with the following input
#RHF/631G(d) stable scf=(tight,qc)
RHF/STO3G 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_{2}, this can immediately be solved by performing an UHF/631G(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/631G(d) guess=(INDO,mix) stable scf=(tight,qc)
UHF/STO3G 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.