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

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.

• "Gaussian" says: Note that analytic frequency calculations are only valid if the wavefunction has no internal instabilities. – pH13 - Yet another Philipp Jun 23 '17 at 15:12
• Yes, this definitely points to its importance. That being said, rarely have I seen people check the stability of the wavefunction, which led me to ask this question since people must have some idea the wavefunction may be unstable in order to check its stability. – Argon Jun 24 '17 at 4:06
• No idea how standard it is, but in my old group, they checked for stability every single time for their open-shell systems, a routine maintenance like checking for imaginary frequencies. These were mononuclear iron or nickel complexes with a ~30 atom, pyridine/phenylene-diamine-based hexadentate ligand. – Blaise Dec 22 '18 at 19:12

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
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)
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
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.
[...]


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.

• Thank you so much! This is extremely helpful and just what I was looking for. I never thought about checking for stability, but all of my systems are open shell, so I suppose I should get in the habit of just checking it seems. – Argon Jun 22 '17 at 18:23
• Speaking anecdotally, I once ran a bunch of calculations on an open-shell nickel system, only to find out after my first 20-30 calculations that half of them weren't stable. The difference post-correction was also quite noticeable. The good thing was that the stable calculations went pretty fast relative to the main calculations, and after correcting an instability the system reliably remained "stable" through further optimization, if that's any consolation. – Blaise Jun 25 '17 at 19:21
• That's exactly what I'm looking to avoid since I'm just starting this project on an open-shell copper system, so the anecdote is quite helpful, thank you for sharing. It will certainly be worth me just running the stability check and playing it safe! – Argon Jun 27 '17 at 23:38