For DFT and Hartree-Fock, how can we know that we have a true minimum? Is there an equivalent to the "frequency analysis" for geometry optimization?

Question

In a geometry optimisation, one can check if they have reached a minimum by checking the phonons. And to a certain extent, the global minimum can be confirmed by relaxing from different starting geometries (sometimes with a bit of initiation).

In SCF calculations (such as Hartree-Fock), is there a way of confirming if one has obtained the global minimum? I assume one could change the initial guess and analyse the converge behaviour, but I was hoping for something a bit more analytical.

Answer 1

You're right about the ability to change the initial guess repeatedly until you get the lowest energy, and this is how it's done in software like MOLPRO which don't offer "stability analysis".

However in software like GAUSSIAN and CFOUR, you can do something called stability analysis, which is described for example in the GAUSSIAN documentation under the command Stable (make sure to click on Options and see all the various different options for using the stability analysis):

"This calculation type method requests that the stability of the Hartree-Fock or DFT wavefunction be tested. Gaussian has the ability to test the stability of a single-determinant wavefunction with respect to relaxing various constraints [Seeger77, Bauernschmitt96](see also [Schlegel91a]). These include:

• Allowing an RHF determinant to become UHF.
• Allowing orbitals to become complex.
• Reducing the symmetry of the orbitals."

In CFOUR you can use the keyword HFSTABILITY=ON to tell CFOUR to do the stability analysis calculation (this basically just prints out the lowest eigenvalues and eigenvectors of the relevant orbital Hessian, but does not "follow" them to try to get a better HF solution), or HFSTABILITY=FOLLOW can be used to actually try to get a lower HF solution if an instability is found during the stability analysis. More about this can be found at the CFOUR documentation for HF stability analysis.

If you properly use the Stable keyword and its relevant options in GAUSSIAN, or the HFSTABILITY=FOLLOW keyword in CFOUR, if lower HF solutions do exist after the first HF calculation is done, then more HF calculations will be done, and the final one should (ideally) give you an energy that can't be beat (within the basis set chosen, as well as the symmetry constraints, although I know CFOUR forces you to use REF=UHF and SYM=OFF so symmetry constraints may not even be possible when doing this).

Answer 2

Just like geometry optimization, there is no practical way to be 100% sure that you have the global minimum of SCF solutions.

But there are checks you can do to make sure that the SCF solution you got is a reasonable minimum. One of them is checking the electronic hessian at the SCF solution and determining the lowest eigenvalues. If one or more negative eigenvalues are found then the solution is a saddle point in the wavefunction space. Read Nike's answer on how to do this with Gaussian and CFOUR.

There are options to do this in other softwares as well, however, not all features might be available.

Orca

In Orca, you can request a stability analysis from %scf block—

! B3LYP def2-SVP NoRI SP 
# only single point calculations supported
%scf
  STABPerform true  # turns on SCF stability analysis
  STABRestartUHFUnstable true # restart the UHF-SCF calculation if unstable
  STABNRoots 3  # (default 3) number of lowest eigenvalues sought
  STABlambda +0.5 # (default +0.5) mixing parameter for new UHF-SCF calculation
end 

*xyz 0 1
H 0.0 0.0 0.0
H 0.0 0.0 1.4
* 

The STABlambda parameter determines how much the old SCF and new orbitals for the new guess. Note that the stability calculation can be done for both RHF/RKS and UHF/UKS references, but restarting the calculation if the solution is unstable (STABRestartUHFUnstable) is only allowed for UHF/UKS references.

Refer to the Orca manual for other available options.

GAMESS

GAMESS can also do stability analysis, but only for RHF references. The keyword UHFCHK=.TRUE. has to be added to $CONTRL section.

! B3LYP/aug-cc-PVDZ on H2
 $BASIS GBASIS=ACCD $END
 $CONTRL SCFTYP=RHF RUNTYP=ENERGY ISPHER=1
 DFTTYP=B3LYP UHFCHK=.TRUE. $END

 $DATA 
Title
C1
H     1.0    -7.56825     2.73003     0.00000
H     1.0    -6.86039     2.74400     0.00000
 $END

This checks pairs of orbitals for RHF->UHF instabilities. By default, it only checks HOMO-1, HOMO, LUMO, LUMO+1 and LUMO+2; but more orbitals can be checked by using the NHOMO and NLUMO keywords in the same section.

You can also turn off any type of symmetry usage by adding the keyword NOSYM=1 in $CONTRL, and see if that changes the SCF solution.

Keep in mind though, that none of these methods are fool-proof and they should not be used as black-boxes. It is always good to check the orbital populations, energies, spin densities etc. to see if that makes chemical sense.

Also note that these codes generally perform calculations with Gaussian type basis sets, and mainly for small molecules. In your question, you used the word "phonon" which makes me think that you might be considering periodic systems with plane wave basis sets. I do not have any experience with those but I believe that the same principles should apply.