So you optimized your molecule, you found all real frequencies, you’re done! right? Well, maybe not. When molecules are well-behaved you probably can either use a restricted or an unrestricted method to find the optimized wavefunction and derive as many properties as you want from it; this occurs for most organic molecules and to some extent to main group elements with little to no radical character. But from time to time, it just so occurs that the wavefunction stemming from your calculation is not the most appropriate to describe your electronic structure. This analysis allows you to quickly check for the wavefunction stability regarding various constraints like allowing a restricted wavefunction to become unrestricted, or allowing orbitals to become complex. In fact, vibrational analysis are only valid if the wavefunction is stable.

Molecules with a biradical character yield instabilities from restricted calculations, also molecules with a triplet nature can show instabilities when calculated at the singlet multiplicity, finally molecules with a multireference character show these instabilities. Therefore, it is important to check the stability of our calculated wavefunctions to fully understand the nature of our molecules.

Let’s take the textbook example of O2, we already know it’s a triplet in the ground state but for the sake of demonstration pretend it’s just another unknown molecule.

I ran a CCSD(T)/aug-cc-pVTZ optimization for the singlet state:

%chk=o2_sing.chk

#p opt ccsd(t)/aug-cc-pvtz freq=noraman

title

0 1
O 
O   1   b1

b1   1.2100

By the way, you need to use a Z-matrix with no more than 50 variables in order for the optimization to carry, otherwise you get the error below, and they need to be variables, it doesn’t work if you place values directly in the Z-matrix. 

 ************************************************
 ** ERROR IN INITNF. NUMBER OF VARIABLES (  0) **
 **   INCORRECT (SHOULD BE BETWEEN 1 AND 50)   **
 ************************************************

From this calculation we get a single real frequency, which could lead us to believe we’re done and ready to calculate further properties. But let’s check the stability of the obtained wavefunction in a second calculation reading it from the checkpoint (we can also make a two step job with the –Link1– option):

%chk=o2_sing.chk

#p geom=allcheck stable guess=read ccsd(t)/aug-cc-pvtz

Be careful NOT to use the stable keyword in your first job, otherwise you’ll analyze the stability of the guess (which is the Harris guess in Gaussian16 by default) and not the final optimized wavefunction.

The output for the wavefunction stability looks as follows:

 ***********************************************************************
 Stability analysis using <AA,BB:AA,BB> singles matrix:
 ***********************************************************************
 1PDM for each excited state written to RWF  633
 Ground to excited state transition densities written to RWF  633

 Eigenvectors of the stability matrix:
 
 Excited state symmetry could not be determined.
 Eigenvector   1:      Triplet-?Sym  Eigenvalue=-0.1434007  <S**2>=2.000
       7 ->  9         0.67311
       7 -> 16        -0.16817
 
 Excited state symmetry could not be determined.
 Eigenvector   2:      Triplet-?Sym  Eigenvalue=-0.1130283  <S**2>=2.000
       8 ->  9         0.67626
       8 -> 16        -0.15201
 
 Excited state symmetry could not be determined.
 Eigenvector   3:      Singlet-?Sym  Eigenvalue= 0.0000259  <S**2>=0.000
       8 ->  9         0.68638
       8 -> 16        -0.13611
 
 Excited state symmetry could not be determined.
 Eigenvector   4:      Triplet-?Sym  Eigenvalue= 0.0778065  <S**2>=2.000
       6 ->  9         0.63124
       6 -> 16        -0.14471
       7 -> 18        -0.21927
 
 Excited state symmetry could not be determined.
 Eigenvector   5:      Singlet-?Sym  Eigenvalue= 0.2344834  <S**2>=0.000
       6 ->  9         0.64723
       6 -> 16        -0.13339
       7 -> 18        -0.19466
 
 Excited state symmetry could not be determined.
 Eigenvector   6:      Singlet-?Sym  Eigenvalue= 0.2612342  <S**2>=0.000
       6 -> 17        -0.10327
       6 -> 18         0.26839
       7 ->  9         0.59483
       7 -> 16        -0.10098
       8 -> 13         0.13453
 The wavefunction has an RHF -> UHF instability.

As we can see, the first section shows a couple of transitions with negative energy to a triplet, which has no physical meaning, but it hints to the possibility of having a triplet wavefunction more stable than the singlet one being analyzed, however further calculations are required to confirm it. The final message reads: “There is a RHF -> UHF instability”, from which, at this point, we can conclude that the restricted version of the employed method is not the best option for the molecule at hand. From here we can take a few directions: changing the multiplicity state to a triplet, or using a Restricted Open shell (RO) method, which keeps as many electrons paired as possible. Keeping the multiplicity state but changing to an unrestricted method makes no difference. So, let’s repeat the process for O2, but lets change the multiplicity state to a triplet from a singlet (for comparison, the following input is divided in two jobs in one file).

%chk=o2_trip.chk
%nprocshared=12

#p opt ccsd(t)/aug-cc-pvtz freq=noraman

title

0 3
O
O 1 b1

b1 1.2100

--Link1--
%chk=o2_trip.chk

#p geom=allcheck guess=read stable ccsd(t)/aug-cc-pvtz

Now, in the output we get the following information in the stability section
 ***********************************************************************
 Stability analysis using <AA,BB:AA,BB> singles matrix:
 ***********************************************************************
 1PDM for each excited state written to RWF  633
 Ground to excited state transition densities written to RWF  633

 Eigenvectors of the stability matrix:
 
 Excited state symmetry could not be determined.
 Eigenvector   1:  3.026-?Sym  Eigenvalue= 0.0219147  <S**2>=2.039
      6B ->  8B       -0.63956
      6B -> 16B        0.27195
      7B ->  9B        0.63956
      7B -> 15B       -0.27195
 
 Excited state symmetry could not be determined.
 Eigenvector   2:  3.026-?Sym  Eigenvalue= 0.0219147  <S**2>=2.039
      6B ->  9B        0.63956
      6B -> 15B       -0.27195
      7B ->  8B        0.63956
      7B -> 16B       -0.27195
 
 Eigenvector   3:  3.063-SGU  Eigenvalue= 0.1798098  <S**2>=2.096
      7A -> 17A       -0.12571
      7A -> 18A        0.30869
      7A -> 32A       -0.11437
      8A -> 12A        0.11528
      8A -> 29A        0.11204
      9A -> 13A        0.11528
      9A -> 30A        0.11204
      5B -> 18B        0.14816
      6B ->  8B        0.56463
      6B -> 16B       -0.22315
      7B ->  9B        0.56463
      7B -> 15B       -0.22315
 The wavefunction is stable under the perturbations considered.

We see that now the resulting wavefunction is stable! Now the energy for the triplet calculation is 55.81 kcal/mol lower than for the singlet state.

Sometimes, using guess=mix is required to find broken symmetry solutions for biradical molecules, because the SCF in Gaussian will try to land a restricted or unrestricted guess from the start. There are several options in the use of the stable keyword but perhaps the most useful is guess=opt which reoptimizes the wavefunction if an instability is found.