From:	"M. Nicklaus" <mn1 %! at !% helix.nih.gov>
Date:	Fri, 1 Nov 1996 19:19:46 -0500 (EST)
Subject:	G:SUMMARY: Wavefunction Instability

Dear CCL'ers,

A few weeks ago I asked the question(s) listed below concerning
wavefunction instability encountered in the calculation of a
cation complexed with the hydroxyl groups of catechol.  I am
summarizing the answers I got, slightly edited for the sake of
brevity.  Thanks to all who responded.

The bottom line seems to be:
1. This instability is probably an artifact of the HF method; and
2. since the energy difference is so small anyway, forget about it.

I performed the UHF/6-311G** geometry optimization in the meantime,
reading in the wavefunction from STABLE=OPT (at the geometry
optimized at RHF/6-31G*).  The energy change this time, i.e. during
geometry optimization, was even smaller (0.25 kcal/mol) than the
change when going from SP --> STABLE=OPT (1.45 kcal/mol).  The spin
contamination increased a bit, from
S**2 before annihilation     0.3882,   after     0.2298
to
S**2 before annihilation     0.4404,   after     0.3042

From these results, and the responses, I tend to deduce that this is
about as good as it gets (at least at the HF level of theory), and
that I can probably reasonably trust the energies I've got so far.
If anyone thinks otherwise, please let me know.

Marc

------------------------------------------------------------------------
 Marc C. Nicklaus                        Lab. of Medicinal Chemistry
 e-mail: mn1 : at : helix.nih.gov               National Cancer Institute, NIH
 Phone:  (301) 402-3111                  Bldg 37, Rm 5B29
 Fax:    (301) 496-5839                  BETHESDA, MD 20892-4255    USA
         WWW:  http://www.nci.nih.gov/intra/lmch/MCNBIO.HTM
------------------------------------------------------------------------

############################ Original Question #########################

I am calculating, at the HF/6-311G** level in Gaussian 94, the divalent cation
Mg++ complexed with the hydroxyl groups of catechol.  Both catechol and the
ion (should) have a spin multiplicity of 1, so that's what I used for the
complex.
   At the current geometry, which was optimized with a smaller basis set, I
performed a wavefunction stability test (keyword "Stable").  I got the message
          The wavefunction has an RHF -> UHF instability.
Running a "Stable=Opt" job produced a stable wavefunction at a slightly lower
energy.  The differences are not huge: the energy is lowered by 1.45 kcal/mol,
and the changes in the charge on the Mg++ ion (both Mulliken and NPA) are less
than 10^-3.  (Interestingly, the complex of Mg++ with the aromatic ring of
the catechol molecule did not show any wavefunction instability.)  Using UHF
instead of (R)HF, while keeping mult. = 1, produced exactly identical results
compared with the HF calculations, both with and without the keyword
"Stable=Opt".
The spin contamination produced by the wavefunction optimization is
 Annihilation of the first spin contaminant:
 S**2 before annihilation     0.3882,   after     0.2298

Before I burn weeks' worth of CPU cycles on potentially unnecessary and/or
meaningless calculations, I'd like to ask for advice on:

(a) Are these results a reason for concern, i.e. does S**2 = 0.3882/0.2298 mean
    that the optimized wavefunction is useless; or, on the other hand, does the
    small energy change of only 1.45 kcal/mol mean that I shouldn't even bother
    optimizing the wavefunction?

(b) Since I intend to optimize the geometry of the complex at the HF/6-311G**
    level, is it possible the the wavefunction instability will go away anyway
    after that?  [I guess I could simply try that, but see (c).]

(c) If I (should) use "Stable=Opt", how do I combine that with the geometry
    optimization?  Can you use both in the same job?  If so, when does G94 do
    the wavefunction optimization; i.e., (only) for the initial geometry,
    and/or for each intermediate geometry, and/or (only) for the final geometry?
    This would appear to have a major impact on the CPU time used.  If you can't
    combine "Stable=Opt" and geometry optimization in a single job, how should
    one proceed?

(d) Any other comments, ideas, suggestions or warnings concerning this problem?

################################# Responses ###############################

From: David Heisterberg 

Do you mean you did a UHF calculation starting from Gaussian's Huckel
guess (or your RHF results)?  And did you use GUESS=MIX?  Anyway, if
you run UHF with the results from the STABLE=OPT run you should recover
the lower energy solution.  I bring that up because what you might want
to do is look at the UHF natural orbitals and I don't know off-hand
whether you can get them as a result of an RHF stability calculation.
The UHF NOs may tell you more about the nature of the instability.
If the interaction between the Mg++ and the OH group is fairly weak at
your current geometry the system may not be well described by a single
determinant.  Look for NOs with populations greater that about 0.02 and
 less than 0.98.

But you may find that the instability goes away as you optimize the
geometry with the larger basis.  You might kill two birds with one
stone by doing a UHF optimization starting from the STABLE output.
If the spin contamination gets small at convergence you're probably
ok.

-------------------------------------------------------------------

From: Frank Jensen  dou.dk>

        Bottom line: don't bother with the UHF wave function.
If you want the singlet, run everything RHF. The UHF instability
is due to a small amount of electron correlation being
included in UHF. If you where to do this at R/U MP2, the
RMP2 would be the lowest energy wave function.

-------------------------------------------------------------------

From: Ulrike Salzner 

HF theory has a tendency to produce triplet (and sometimes singlet)
instabilities that are artifacts of the method. That can be
understood by the neglect of part of the electron correlation in HF
theory. Since the Pauli principle is obeyed by using determinatal
wave functions, electrons with same spins are correlated. They can
not be at the same place. Electrons with different spins, however,
are not correlated. That is, the error in the HF approximation is
smaller for triplets than for singlets. This can lead to erronous
preference for the triplet.

Since your system should be a singlet, as you said, and since the
energy differences between your singlet and triplet calculations are
quite small, I think you can forget about the triplet instability.

-------------------------------------------------------------------

From: Doug Fox 

    I suspect that the reason you got back the same energy is because
you did not read in the result of the STABILITY calculation as the initial
guess.  The default initial guess is RHF so just saying UHF is not
enough to break the spin symmetry.  You might try UHF with GUESS=MIX
or read in the solution from STABLE=OPT.

    Will this instability go away, no option but to check.  Normally
I would not expect so unless the structure changes dramatically and
what you are seeing is due to a bad initial guess at the structure.
Since many people do optimizations without ever checking it is not
easy to guess how common that is.  So I would optimize normally
and then run STABLE=OPT again at the end to see if you have stayed
on the more stable surface.

-------------------------------------------------------------------

From: Michael Braeuer 

I had to deal with the same problem optimizing a simple problem as acrolein.
So far I didn't find any explanation but I'm working on it.

-------------------------------------------------------------------

From: Anthony P Scott 

I would be interested in the general comments on this. I have noticed
similar spin contamination values and have not been to sure what to make
of them

I think that an instability in the wavefunction (given a suitable basis)
will not go away with increasing basis set size. I would be happy to be
corected on this however.,

I would proceed (and do) as follows.
a. optimize the geometry at rhf or whatever.
b. do a stability run with stable=opt.
c. with this wavefunction (guess=read) reoptimize the geometry.

Generally you will find that the geometry does not change much.

######################## End of Responses #######################


Similar Messages
10/01/1996:  G:Wavefunction Instability
08/01/1995:  Spin contamination, effect on energy and structure.
12/16/1994:  Spin contamination & AM1 "ROHF" versus UHF
04/23/1992:   Huckel MO Theory software
08/29/1995:  Summary for "Calculation on biradical"
08/01/1996:  Re: CCL:M:Heat of formation calculation using MOPAC.
07/22/1992:  Re: gaussian wavefunctions
12/09/1997:  CCL:G:Summary - open shell calculation
03/11/1996:  Law of conservation of difficulty: violations.
02/01/1997:  G94(Sum.+Q): Problem with PES Scan in DFT


Raw Message Text