When is it necessary to check wavefunction stability in density functional theory?
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 -&gt; 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.