So, you’ve optimized your molecule, obtained only real frequencies, and everything looks perfect—done, right? Well… not necessarily.
For many well-behaved systems—especially typical organic molecules or
main-group compounds with little to no radical character—a standard restricted
or unrestricted calculation is sufficient. In such cases, the optimized
wavefunction is usually reliable, and you can confidently extract molecular
properties from it.
However, there are situations where the computed wavefunction is not
the most appropriate representation of the true electronic structure. This
is where wavefunction stability analysis becomes essential. It allows
you to test whether your solution remains stable when certain constraints are
relaxed—for example, letting a restricted wavefunction become unrestricted or
allowing orbitals to become complex. Importantly, even vibrational analyses are
only strictly valid if the underlying wavefunction is stable.
Instabilities are common in systems with:
- Biradical character (often problematic for
restricted methods)
- Incorrect spin states (e.g., singlet calculations for
inherently triplet systems)
- Multireference character, where a single-determinant
description is insufficient
For such systems, checking wavefunction stability is not optional—it is
crucial for correctly describing the molecule.
A Classic Example: O₂
Let’s consider molecular oxygen as a textbook example. We already know
that its ground state is a triplet, but suppose we treat it as an unknown
system and perform a singlet calculation.
After optimizing the singlet structure (e.g., at the CCSD(T)/aug-cc-pVTZ
level) and confirming that all frequencies are real, one might be tempted to
proceed. However, before doing so, we should test the stability of the obtained
wavefunction.
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
This is done in a second calculation using:
%chk=o2_sing.chk#p geom=allcheck stable guess=read ccsd(t)/aug-cc-pvtz
A key point: do not include the stable keyword in the initial job,
as this would analyze only the initial guess (typically the Harris guess in
Gaussian), 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.
Interpreting Stability Analysis
In the stability output, negative eigenvalues indicate instabilities. For
the singlet O₂ calculation, you may observe negative eigenvalues corresponding
to triplet excitations, followed by a message such as:
“The wavefunction has an RHF → UHF instability.”
This tells us that the restricted singlet solution is not stable and that
a lower-energy unrestricted solution likely exists. In other words, the
electronic structure is not properly described at this level.
At this stage, you have a few options:
- Change the spin multiplicity
(e.g., singlet → triplet)
- Use a restricted open-shell
(RO) method
- Explore unrestricted solutions
more carefully
Simply switching from restricted to unrestricted without changing
multiplicity will not resolve the fundamental issue.
Correcting the Problem
If we repeat the calculation for O₂ with triplet multiplicity, and
then perform the same stability analysis, we now obtain:
“The wavefunction is stable under the perturbations considered.”
This confirms that the triplet solution is the correct and physically
meaningful description. Notably, the triplet state is also significantly lower
in energy (by ~55.8 kcal/mol) compared to the singlet, reinforcing this
conclusion.
%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.
Practical Tips
- For biradical systems, you
may need to use:
guess=mix
to obtain broken-symmetry solutions, as Gaussian often defaults to
restricted-like guesses.
- A very useful option is:
stable=opt
which automatically re-optimizes the wavefunction if an instability is
detected.
Take-Home Message
Obtaining an optimized structure with all real frequencies does not
guarantee that your job is complete. A stable geometry built on an unstable
wavefunction can lead to incorrect interpretations and unreliable properties.
Wavefunction stability analysis is therefore a critical final
checkpoint—ensuring that your electronic structure is not only mathematically
converged, but also physically meaningful.
Comments
Post a Comment