When is it necessary to check wavefunction stability in density functional theory?

Question

In Gaussian, there is a stable keyword that checks the stability of the wavefunction. Using stable=opt reoptimizes the wavefunction until a stable one is found if there is an instability.
My question is - when do you use this stability check? Presumably one must have some hunch in advance that the wavefunction may be unstable in order to use this keyword, but what insight would lead one to this conclusion?

To add some context, this question came to mind after following instructions on how to model a system with antiferromagnetic coupling (e.g. as shown here). When I tried this procedure on cupric acetate hydrate, which has a dinuclear structure of copper(II) atoms and is known to be antiferromagnetic, I had to use the stable=opt keyword to stabilize the wavefunction resulting from the "fragment guess job" otherwise it erroneously appeared that the antiferromagnetic singlet was higher in energy than the triplet. This made me wonder when I should be using stable=opt in general (perhaps outside the context of antiferromagnetic coupling calculations) and also what it is specifically doing.

Answer 1

The very short answer is: whenever you are in doubt.

The stability check is probably unnecessary when you are dealing with closed shell molecules, where you can be pretty sure that a restricted methodology is reasonable. For example, when you calculate water at the RHF level of theory, you won't really need to look for a more stable wave function, because it would not be physically sensible. An UHF approach in principle will still lead to a lower energy since it has more flexibility.

Whenever you are calculating molecules in which you would assume some open shell character you should be considering checking wave function stability. It usually doesn't hurt to perform one check too many. Discovering instabilities later can be a lot more frustrating. That does not necessary mean that you are gunning for the correct ground state, as you are still restricted on spin. I'll include a small, but well-known example below: $\ce{O2}$. I am using Gaussian 09 Rev. D.01 for the calculation, but it should work with anything more recent. (It won't work the same way with anything older.) Also note, that I am not conducting geometry optimisations, but use $\mathbf{d}(\ce{O-O})=\pu{120pm}$ to keep it simple.

First let's run an ordinary, restricted calculation. If you are not using the SCF(xqc) keyword, you should immediately notice, that something is wrong.

%chk=bp86svp.rhf.sing.chk
#p RBP86/Def2SVP/W06
DenFit             ! Use density fitting to spead up calculation
gfinput            ! Prints basis set in a form suitable for use as general basis set input
gfoldprint         ! Prints the basis set information in the Gaussian format
iop(6/7=3)         ! Switches on printing of all MO
symmetry(loose)    ! Use correct symmetry
scf(xqc,conver=10) ! Quadratic convergence after conventional (for tough cases)
INT(SuperFineGrid) ! Largest grid, most accurate

title card unused

0 1
O  0.6 0.0 0.0
O -0.6 0.0 0.0

You'll find a summary of the energies at the end, if you want to check and run your calculations.

Next up, we'll check and optimise the wave function.

%chk=bp86svp.rhf.sing.stab.opt.chk
%oldchk=bp86svp.rhf.sing.chk
#p RBP86/Def2SVP/W06
DenFit
geom=check guess=read
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
scf(xqc,conver=10) INT(SuperFineGrid)
stable=opt

title card unused

0 1

You should find the following in your output:

 The wavefunction has an RHF -> UHF instability.

Followed by another SCF run. Pay close attention to the following:

 SCF Done:  E(UB-P86) =  -150.201660211     a.u. after    6 cycles
            Convg  =    0.3033D-11                    52 Fock formations.
              S**2 =  1.0016                  -V/T =  2.0049
 = 0.0000 = 0.0000 = 0.0000 = 1.0016 S= 0.6187
 = 0.000000000000E+00
 Annihilation of the first spin contaminant:
 S**2 before annihilation     1.0016,   after     0.0126

You are still finding a singlet state, but now as a UHF wave function. You should later find:

 The wavefunction is stable under the perturbations considered.

What you have found is a broken symmetry singlet. We do know that the ground state should be a triplet. Therefore we'll perform a restricted open shell calculation. We already know this one will have and RHF -> UHF instability, and since we cannot use the quadratic convergence, we cannot use stable either. But it'll provide us with a good starting point for the UHF calculation, or at least an energy to compare others to.

%chk=bp86svp.rohf.trip.chk
%oldchk=bp86svp.rhf.sing.chk
#p ROBP86/Def2SVP/W06
DenFit
geom=check guess=alpha
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
scf(conver=10) ! Quadratic convergence not available for RO
INT(SuperFineGrid)

title card unused

0 3

The important part here is, that it generates a pure spin state wave sunction:

 = 0.0000 = 0.0000 = 1.0000 = 2.0000 S= 1.0000

Another calculation is using the the mixing of HOMO and LUMO to generate a broken symmetry ansatz wave function directly. You'll find the same result as with the approach above.

%chk=bp86svp.uhf.sing.chk
%oldchk=bp86svp.rhf.sing.chk
#p UBP86/Def2SVP/W06
DenFit
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
geom=check guess=(read,mix)
scf(xqc,conver=10)
INT(SuperFineGrid)
stable=opt

title card unused

0 1

You'll therefore find again:

 SCF Done:  E(UB-P86) =  -150.201660211     A.U. after   11 cycles
            NFock= 11  Conv=0.88D-12     -V/T= 2.0049
 = 0.0000 = 0.0000 = 0.0000 = 1.0016 S= 0.6187
 = 0.000000000000E+00
 KE= 1.494731832724D+02 PE=-4.114759744301D+02 EE= 8.357834648903D+01
 Annihilation of the first spin contaminant:
 S**2 before annihilation     1.0016,   after     0.0126

The wavefunction is stable under the perturbations considered.

Now for the last part, run a real UHF calculation:

%chk=bp86svp.uhf.trip.stab.opt.chk
%nproc=8
%mem=8000MB
%oldchk=bp86svp.rhf.sing.chk
#p UBP86/Def2SVP/W06
DenFit
geom=check guess=read
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
scf(xqc,conver=10)
INT(SuperFineGrid)
stable=opt

title card unused

0 3

You'll find as expected

 Annihilation of the first spin contaminant:
 S**2 before annihilation     2.0031,   after     2.0000

 The wavefunction is stable under the perturbations considered.
 [...]
 The wavefunction is already stable.

For giggles, why not consider a quintet?

Command file              Functional            Energy / Hartree ( cycles )
bp86svp.rhf.sing          E(RB-P86)       =       -150.155030787 (      4 )
bp86svp.rhf.sing.stab.opt E(UB-P86)       =       -150.201660211 (      6 )
bp86svp.rohf.trip         E(ROB-P86)      =       -150.213771982 (     10 )
bp86svp.uhf.sing          E(UB-P86)       =       -150.201660211 (     11 )
bp86svp.uhf.trip.stab.opt E(UB-P86)       =       -150.215679340 (      9 )
bp86svp.uhf.quin.stable   E(UB-P86)       =       -149.651916836 (      9 )

In conclusion, whenever you suspect some open shell character, check for stability. If you have well-behaved molecules (singlet), like most organics are, you can probably skip it.
Be aware that a stability optimisation will not necessarily give you the ground state.