Subject: Summary: wavefunction instability vs. spin contamination

From: Theo de Bruin <[email-masked]>
Date: Mon, 18 Sep 2000 12:15:29 +0000

Dear CCL'ers

Some time ago I posted a question about a wavefunction instabiltiy versus
spin contamination problem.  Although I received some replies, I was not
able to solve the problem: i.e. to find a stable wavefunction that yields
good geometries (with respect to crystal data structures).
I can add that the use of a larger basis set does not make the instability
disappear (as was suggested by one the repliers). Apart from using
d-functions on C and O atoms, we now also added f-functions on the Co atoms
(these atoms contain the highest spin density using UB3LYP) but this had no
influence on the (in)stability of the wavefunction.
Anyway, I would like to thank Joe Harrison, Robert Szilagyi and John
McKelvey for their help.

Here is the orginal question, followed by the replies:

*********
QUESTION
*********
At the moment we are performing calculations on several (relatively large)
di-cobalt complexes using the B3LYP functional and LanL2DZ basis set.
Restricted B3LYP calculations yield geometries that fit quite well the
corresponding crystal structures.  However, from computations with the G98
keyword STABLE it was clear that wavefunctions showed a RHF --> UHF
instability. For example:  Eigenvector   1:   Triplet-?Sym  Eigenvalue =
-0.0067007  (it seems that the negative eigenvalue is rather small....?)
Performing the same calculations with the UNrestricted B3LYP (using
guess=mix and changing the DIIS for the QC algorithm) does yield stable
wavefunctions, which are about 0.1 to 0.2 kcal/mol lower in energy.
However, the obtained geometries with the unrestricted wavefunctions
deviate substantially from the crystal structures (delta r(Co-Co) > 0.2
angstrom!), and moreover there is a large spin contamination (S**2 is up to
0.74 for a singlet ground state), which probably results from the (too)
large Co-Co distance.

We have calculated that the triplet electron configuration has a much
higher energy (> 50 kcal/mol).  Unfortunately, we were not able to perform
successful CASSCF calculations to find out if there was contribution of the
triplet in the ground state, since the necessary number of electrons and
orbitals in the active space was far too large.

Since MP2 calculations on 1st row transition metals are unreliable, higher
(than HF or DFT) correlations methods are completely impractical for this
size of systems, and the instability seems to be very small (small neg.
eigenvalue and small energy difference between restr. and unrestr.
wavefunctions), it is very tempting to continue the work with the
restricted B3LYP method.  Especially, because we are mostly interested in
the geometries and their relative energies.

Does anybody have some tips or advice how to continue?

********
ANSWERS
********
Joe Harrison wrote:
I'm not sure if you are using a package that would allow you to look
at the spin densities graphically, but if so, I would recommend looking
at the difference in spin densities for outer atoms of your complex.
I suspect that your UHF result allows for "surface" spin density to
build up whereas your RHF result suppresses that.  That would seem to be
consistent with the increase in Co-Co distance (movement of the electron
density from between the atoms to the "surface").  If your objective is
to use your RHF calculation as a paradigm for simulating a crystal, then
>from that perspective, the use of UHF would be analogous to allowing
surface states to arise.  If this were a cluster simulation of a covalent
crystal, the effect you describe would be like the problem of "dangling
bonds" in those applications.  Perhaps there would be something
analogous to hydrogen terminators that could be used to suppress the
tendency for the spin density to migrate to the outside of the complex.
However, hydrogen terminators can result in a monumental explosion of
the size of the calculation, so whatever the analog might be, it would
probably suffer the same drawback.  If nothing else, if what I have
suggested is in fact what you find, you could use it as partial justi-
fication for the continued use of RHF.  BTW if you don't have the
ability to look at the spin density, the spin populations in the
G98 output might also give such information with the usual caveat
in regard to basis sets with diffuse orbitals.

Robert K Szilagyi wrote:
my guess would be carefully check the wavefunctions. You can easily
analyse the gross population analyses, spin density of each orbitals.
Also you might consider different spin-coupling schemes (ferromagnetic,
antiferromagnetic) and try to calculate the broken symmetry wavefunction
where you mix the singlet and triplet wavefunctions. I don't know your
system, whether it is feasible or not!
I strongly recommend to correlate these data with experimental values of
excitation energies, ground state spin densities, etc.
Once you have the right wavefunction, the geometry of the complex should
be okay as well and not other way around!

John McKelvey wrote:
While Post-Docing with Prof Berthier someww time ago (> 20 years) I
discovered that UHF triplet spin contamination of singlet molecules seems
to disappear as the quality of the basis set is improved.  Later someone
else reported that the
contamination decreases even further if the geometry is simultaneously
optimized...

-----------------------------------------------------------------------
Dr. Theodorus J.M. de Bruin       Phone  +33 4 76 63 56 28
Université Joseph Fourier         Fax    +33 4 76 51 43 82
LEDSS VII  Chimie Théorique       E-mail [email-masked]
UMR CNRS 5616
Bâtiment 52 - Chimie Recherche
DU BP 53
38041 Grenoble Cedex 9
FRANCE
-----------------------------------------------------------------------