Subject: Summary of 'spin contamination'

From: "Michael D. Bartberger" <[email-masked]>
Date: Tue, 1 Aug 1995 14:07:34 -0400 (EDT)

Hello all:

A while back I posted a question regarding the reliability of energies / 
properties computed with UHF where the S^2 value was fairly high ("spin 
contaminated") but the value of S^2 after projection of the first spin 
contaminant was fairly tolerable.  I received quite a few responses and 
would like to thank those who took the time and had the interest to 
respond.  What follows is a summary of responses.  In addition, there was 
a thread a while back on the use of ROHF vs. UHF, which I have also 
included in this post.  The whole body of the message shouldn't be much 
more than 44K or so.

Since the posting of this message, I have come across a paper in J. Phys 
Chem (1995, vol. 99, 8582) by Wong and Radom which is a quite thorough 
theoretical treatment of some radical additions to alkenes, within which 
problems regarding spin contamination are addressed.  This might be of 
interest to those who responded to my original post.

Again, many thanks to the respondees.  

-Michael
________________________________________________________________________________

 Michael D. Bartberger         [email-masked]        TEL: (904) 392-3580
 Department of Chemistry       [email-masked]         FAX: (904) 846-0296
 University of Florida
 Gainesville, FL  32611
 USA
________________________________________________________________________________

******************************

***********************

From [email-masked] Aug  1 13:42:56 1995
Date: Wed, 5 Jul 95 13:31:22 EDT
From: [email-masked]
To: [email-masked]
Cc: [email-masked]
Subject: RE: CCL:Spin contamination

        If memory does not fail me, projection out of the first spin con-
taminant leaves you with an improved set of density matrices for further use
in properties calc-ns.  However, these matrices no longer correspond to the
energy that was determined in the iterative step.  The discrepancy may be
quite significant.  For a *fixed* geometry, one can recalculate the energy
with the improved density matices.

        If geometry optimizations are involved, one has to be concerned that
the geometry obtained without spin-projection may differ in non-trivial
ways from that which would be obtained using spin projection.  There may be
some studies on this, but I cannot call them to mind.

        I think you can see that by the time one gets to energy *differences*,
we have the potential for a mess: mismatches in spin contamination, structures
that have differing energy relationships to their (spin-projected) minima, etc.

        Projected-UHF calculations we in vouge at one time, but I do not
believe that these are incorporated into G92 or G94.  

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                           FREDERIC A. VAN-CATLEDGE

Scientific Computing Division         ||   Office: (302) 695-1187 or 529-2076
Central Research & Development Dept.  ||          
The DuPont Company                    ||      FAX: (302) 695-9658
P. O. Box 80320                       ||
Wilmington DE 19880-0320              || Internet: [email-masked] 
--------------------------------------------------------------------------------
Opinions expressed in this electronic message should ***> NOT <*** be taken to 
represent the official position(s) of the DuPont Company.  

*****> ANY OPINIONS EXPRESSED ARE THE PERSONAL VIEWS OF THE AUTHOR ONLY. <*****
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

************************

From [email-masked] Aug  1 13:43:33 1995
Date: Wed, 5 Jul 95 15:06:00 EDT
From: "Doug S. Dudis MLBP CS" <[email-masked]>
To: [email-masked]
Subject: rohf hessian

Michael,
I'm 99% sure that GAMESS (Iowa State) DOES have analytical hessians
for ROHF hessians.  We use it more than we use G92, since lots of our
problems would have large spin contaminations (yours seem small by
comparison) and must use ROHF.  I agree with your concerns.
Hope it helps
Doug
Polymer Branch
Wright Patterson AFB

*****************************

From [email-masked] Aug  1 13:44:20 1995
Date: Thu, 6 Jul 1995 13:11:56 +0200 (DFT)
From: Mathias Brock <[email-masked]>
To: [email-masked]
Cc: Mathias Brock <[email-masked]>
Subject: Re: Spin contamination

Dear Michael,

this is not an answer to your question (what's behind the G92-
spin annihilation algorithmus?) , but the trouble is not new to me.
Thus, I would be interested in the answers, too.

As far as I understood electron correlation was taken into account at
UHF levels so to vary the spatial coefficient of paired alpha and
beta electrons. 
In a case of high spin contamination on aromatic molecule, I got
a useless structure UMP2; which became only reasonable using cas
method, which was suitable for this PI-system. This example has
shown, that the order of alpha and beta MO's is not the same for
the frontier orbitals,if the difference in the spatial extent of
corresponding electrons is high. 
So just a suggestion: the extent of (a,b)-spatial difference may
by limited, so a useful amount of uhf-corelation can be saved
without unacceptable spin**2 value.

Please send my  the list of answers, if not send to CCL

Mathias

*******************************

From [email-masked] Aug  1 13:44:28 1995
Date: Thu, 6 Jul 1995 16:36:32 +0200 (MET DST)
From: Hans Ueli Suter <[email-masked]>
To: [email-masked]
Subject: Spin-Contamination

Dear Michael,

indeed you have a large spin-contamination. How large the effect is on
your geometries and frequencies is hard to estimate, however.
The annihilation procedure will not help in your case, because, I think
the gradients and hessians are calculated with the not-annihilated
wave-function (better check the manual, probably in Gaussian99...) and
I guess you will not like to do the hessian by hand with the plotted
PUHF-Energy. The advice I would give you, yes it is not very complicated,
use better wave-functions for your molecules. Since you say, they are to
large to do the only reasonable thing (MCSCF-CI!) you should use density
functional theory, probably the Becke3LYP functional. The amount of
computing time (using gaussian) should not be more than 3 times your
UHF calculation. Normaly, there the spincontamination is smaller.

                                           with best regards,

                                              Hans Ulrich  Suter

****************************

From [email-masked] Aug  1 13:45:01 1995
Date: Thu, 6 Jul 1995 17:52:57 +0100
From: Thomas Bally <[email-masked]>
To: "Michael D. Bartberger" <[email-masked]>
Subject: Re: CCL:Spin contamination

Dear Michael,

  annihilation of the first spin contaminant (i.e. the quartet contribution
in your case of a doublet) has no repercussion on the geometry. All it says
if it makes <S**2> drop to an acceptable value is, that there are no impor-
tant contributions to higher spin contaminants (hextets, octets ...).

  If your <S**2> is 1.0 for your radicals, then I would be sceptical with
regard to your geometries because the UHF model will inherently favor
structures with high spin contamination (those lower the energy!). What
is worse, it makes MP2 calculations meaningless for all practical pur-
poses (this is where spin projection comes in: if your projected wave-
function has <S**2>=0.76 then your PUMP energy is probably O.K., but I
think you will not do MP2 on C6F11).

 But why not use ROHF? Gaussian is not very strong in this department but
what keeps you from using GAMESS? It's free and has ROHF analytical second
derivatives plus very efficient ROMP2 code, albeit only for single points.
Its only serious disadvantage in my view is, that it will not let you use
Gaussian-style Z-matrices unless you use no dummies which is nearly impos-
sible if you have rings and give up on keeping control over certain inter-
nal coordinates (well, lest I get flamed by Mike Schmidt, I should say that
*in principle* it is possible, but in practice it is very tedious to say the
least). But otherwise it is a fine program, we use it a lot.

 If you think spin polarization is important for you then you should know
that UHF usually overestimates this something awful (even more if higher
spin contaminants are important). A relatively inexpensive "safe" way to
account for spin polarization is by MCSCF. However, you will want to look
for some guidance there because if you compare different points on a poten-
tial energy surface, then your choice of active space becomes very critical
if you want to avoid discontinuities on your reaction path. This is not
your "black box" type theory!

 However, you are a few steps from one of the most famous schools of quantum
chemistry (Quantum Theory Project) and I am sure Mike Zerner, Rod Bartlett or
some other practically minded person there would be able and willing to
help you with this, if you elect to go this way. Give them my regards if you
do!

 best wishes

  thomas

*-------------------------------------------------------------------------*
|  Prof. Thomas Bally                 |  E-mail:  [email-masked]   |
|  Institute for Physical Chemistry   |  WWW page:                        |
|  University of Fribourg             |  http://sgich1.unifr.ch/pc.html   |
|  Perolles                           |                                   |
|  CH-1700 FRIBOURG                   |  Tel:     011-41-37 29 8705       |
|  Switzerland                        |  FAX:     011-41-37 29 9737       |
*-------------------------------------------------------------------------*

**************************

From [email-masked] Aug  1 13:45:20 1995
Date: Fri, 7 Jul 1995 20:49:43 -0500 (CDT)
From: [email-masked]
To: "Michael D. Bartberger" <[email-masked]>
Subject: Re: CCL:Spin contamination

Michael,
> 
> In G92, after the value of <S^2> is printed in a UHF calculation, there 
> is a value printed after this which shows <S^2> "after annihilation of 
> the first spin contaminant." 
> 
> Could someone kindly explain the significance of this in terms of the 
> "usefulness" of the energy or geometry obtained with the 
> spin-contaminated wavefunction?
> 
   The projected value of <S**2> tells you what happens after projecting the
S+1 state (in your case the quartet). The fact that it gets close to 0.75
would tend to indicate that higher spin states are of limited importance.
However, your geometry, energy, and population analysis are all for the UHF
wavefunction. If you want the spin projected energy, you need to do an MP2
calculation, which in Gaussian will also print the PUHF energy. If you would
like to see the population analysis of the projected (annihilated is actually
more accurate in this case, since it is only annihilating the S+1 state --
see below) you can specify iop(6/15=2) in the route card, telling the
population analysis link to use the projected density. If you want to
optimized the geometry at the projected level, you will have to do it by hand
(and it won't be fun . . . ). Note that spin projection gives you a more pure
spin wavefunction, but one that is necessarily constructed from the UHF
orbitals (that are obviously not very good!). PMPn energies can be
calculated by doing MP4 calculations (see Schlegel paper in JPC, 1988 I
believe -- sorry, I'm at home, so reference not in front of me.)

   Finally, spin contamination is an artifact of the UHF formalism trying
to build electron correlation into the wavefunction by separating the spatial
components of the alpha and beta orbitals. MCSCF methods would be a good
alternative, since symmetry is easier to maintain and correlation
(non-dynamic, anyway) is handled correctly.

Hope this helps,

Chris

-- 

Christopher J. Cramer
University of Minnesota
Department of Chemistry
207 Pleasant St. SE
Minneapolis, MN 55455-0431
--------------------------
Phone:  (612) 624-0859 || FAX:  (612) 626-7541
[email-masked]
http://dionysus.chem.umn.edu/

*************************************************************************
*                                                                       *    
*   The following posts were not in direct response to my question,     *
*   but rather part of a discussion of "ROHF vs. UHF" which transpired  *
*   a while back.  I thought this would be of interest to those who     *
*   were also interested in my question and was worth posting.          *
*                                    -MDB                               *
*************************************************************************

>From [email-masked]  Fri Dec 16 14:22:09 1994
Received: from ohstpw.mps.ohio-state.edu  for [email-masked]
	by www.ccl.net (8.6.9/930601.1506) id OAA25185; Fri, 16 Dec 
1994 14:10:19 -0500 From: <[email-masked]>
Date: Fri, 16 Dec 1994 14:10:00 -0500 (EST)
Subject: Spin contamination & AM1 "ROHF" versus UHF
To: [email-masked]

Dear Netters,

I posted a number of questions to the CCL about spin contamination 
on the 14th November.  Below is the original questions and a summary
of the responses.  I also posted on the 7th of December questions 
relating to the 5th question below.  The responses to this posting 
are also summarized here.

======================================================================
                         14th November Posting
======================================================================

I have a number of questions regarding the effects of spin
contamination on total electronic energies and structures.  My
understanding of spin contamination is that unrestricted Hartree-Fock
(UHF) wave functions are not eigenfunctions of the total spin
operator, so the electronic wave function of interest can be
contaminated by functions corresponding to states of higher spin
multiplicity.  This brings me to several of my questions:

1.  Given that, (a) we have performed an ab initio study on an open
shell molecule using a single-determinant wave function (i.e.,
variational), (b) we have found its theoretical equilibrium geometry
for the ground electronic state, and (c) have projected out ALL
contaminating higher spin multiplicity states along the way to the
minimum energy geometric configuration.  Will the total calculated
electronic energy be the lowest possible calculated energy for the
given basis set and Hamiltonian?

2.  What can be said about 1, above, regarding the total calculated
electronic energy when a perturbation treatment to the configuration
interaction is introduced, e.g., MP4(SDTQ)?

3.  What can be said regarding the total calculated electronic energy,
if in the case of 2, above, condition (1c) is not fully met, i.e.,
some spin contamination remains while going to the optimum geometry.

4.  Given that different spin states correspond to different
equilibrium geometries, is an optimized structure, which had
significant spin contamination in its electronic wave function all the
way down to its "minimum" energy geometrical configuration, some kind
of mixture of structures involving the state of interest and the
different contaminating higher spin multiplicity states?

5.  Regarding semiempirical calculations, in the paper of Novoa et
al., Inorganica Chemica Acta, Vol. 198-200 (1992) 133, the heats of
formation of some very large carbon clusters were calculated.  In
their paper the authors state:  "For an odd-membered linear C_n with
13 <= n < 20, the AM1 calculations predict the triplet state to be
more stable than the singlet, due mainly to the spin contamination of
the UHF calculations."  The calculated stabilities of these larger
linear carbon clusters are not expected based on what is known for the
smaller linear clusters where, for odd n (n > 1), the singlet states
are more stable than the triplet states.  My question is this, is it
possible that the authors are not correct regarding their reason for
the for greater stabilities of the triplet states?  (It is not my
intention to attack the above authors work.  I merely wish to evaluate
the quoted heats of formation because I require some kind of estimate
for the heats of formation for these species.)

Regards,

Ryan Bettens,
OSU Physics Department,
[email-masked]

                       ________________________
----------------------| Responses to Questions |----------------------
                       ------------------------

From:  Christopher J. Cramer,
       University of Minnesota,              
       Department of Chemistry               

   Re your queries on spin contamination. I can not speak to all of 
your questions, but one point I think worth making is that the UHF SCF
procedure is variational in that it produces the unique orthogonal set
of MO's that give rise to the lowest energy wavefunction for whatever 
basis set has been chosen. However, those MO's derive from the UHF 
calculation itself, which includes all spin contamination. Any 
projection operator formalism applied thereafter can provide you some 
(significant) improvement in energy because of the annihilation of 
higher spin states, BUT that procedure does NOT reoptimize the MO's. 
So, even if you optimize a geometry by hand using PUHF (or PMPn) 
energies (since gradients for these methods are not, to my knowledge, 
available) you will still suffer from potentially poor wavefunctions.

   Alternatives include ROHF (which I personally don't like because we
tend to be interested in spin-polarization, which you can't get at 
this level), spin-adapted MCSCF treatments, and, something you might 
consider, spin-polarized DFT. The latter is in essence a UHF-like 
calculation using an approximate scheme for exchange/correlation (as 
all DFT does), but it has been shown to be far less prone to spin 
contamination.

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

From:  Dave Ewing,
       John Carroll University,              

See warren J. Hehre, et. al., Ab initio Molecular Orbital Theory (John
wiley & Sons, 1986), pp. 203-4.

My own experience has been that structures are OK as long as the spin
contamination isn't too gross, e.g. no more than <S2>=1.0 for a 
doublet.

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

From:  I. Shavitt,
       Ohio State University,              
       Department of Chemistry               

1.  The wave function obtained by projecting out all spin
    contamination AFTER having calculated a UHF solution
    is not the lowest-energy solution for the spin-state
    in question.  This is so because the orbitals have been
    optimized for the spin-contaminated function, before
    projection, not for the projected function.

2.  After the addition of a correlation treatment, like
    MP4(SDTQ), the effects of spin contamination would
    usually be diminished, though not eliminated.  The
    effect on the energy cannot be predicted with confidence,
    because there are diverse factors determining how the
    choice of orbitals will affect the correlation energy.
    In fact, orbitals which give the lowest Hartree-Fock
    energy are not guaranteed to produce the lowest
    correlated energy, though they usually do.  Coupled
    cluster methods, such as CCSDT, are less sensitive to
    the choice of orbitals, and therefore less likely to
    suffer from the use of nonoptimal orbitals.

3.  See 2.

4.  The effect of spin contamination on optimized geometries
    depends both on the degree of spin contamination and on
    the relative energies and characters of the higher-spin
    states.  I don't think I can give a general answer to this
    question.

5.  Personally, I am very skeptical about the ability of
    methods like AM1 to give definitive answers to questions
    concerning the relative energy of different electronic
    states.  But the same is true, even more strongly, for
    UHF calculations.  For some molecules it is very difficult
    to determine which structure or multiplicity is lower in
    energy.  An example is C_4, for which a linear triplet state
    and a cyclic singlet are almost of the same energy, and
    different levels of even high-quality calculations give
    conflicting answers.  My impression was that odd-membered
    carbon chains have a singlet ground state, and it would take
    much more than AM1 calculations to convince me otherwise.

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

From:  Ab Initio Molecular Orbital Calculations for Chemists, 2nd ed.,
       W.G. Richards and D.L. Cooper,              
       Clarendon Press, Oxford, 1983.               

I found this comment in the above book after posting the above 
questions.  The comment essentially answers my question 4.  I quote 
this book here (from page 59) as it may be of convenience to people.

     The average interaction interaction for alpha-spin electrons
     and beta-spin electrons may be different in open-shell 
     systems so that it is not unreasonable that orbitals 
     differing only in their spin quantum number should have 
     different spatial functions.  This is the origin of the
     unrestricted Hartree-Fock (UHF) method implemented in some 
     programs.  Unfortunately, UHF wavefunctions are not always 
     eigenfunctions of S^2.  If for example we carried out 
     calculations on a molecule with a doublet state and a 
     quartet state that were close in energy, then the 
     equilibrium geometry of the doublet state would be partly
     characteristic of that for the quartet state.
     
     Although there are methods for projecting out spin 
     eigenfunctions at the end of the UHF calculation, many 
     theoreticians prefer the restricted Hartree-Fock (RHF) 
     method and this is implemented in many open-shell SCF 
     programs.  In the RHF method, orbitals differing only in 
     spin quantum number . . . have identical spatial parts.

======================================================================
                         7th December Posting
======================================================================

Thanks to those of you who responded to my last posting on the 14th of
November regarding spin contamination.  Before summarizing these
responses I would like to ask one further question which is related to
my last question in the previous posting.  I will then compile all
responses and post a summary.  My question concerns ROHF versus UHF
calculations using the AM1 semiempirical model.  Is anybody aware of
any published work which recommends the use of either in calculating
heats of formation, specifically for hydrocarbon open shell species?

I have examined the original paper of Dewar, M.J.S.; Zoebisch, E.G.;
Healy, E.F. and Stewart, J.J.P., J. Am. Chem. Soc., 107 (1985) 3902
where AM1 was introduced.  In this work the heats of formation of 6
hydrocarbon radicals (all doublets) where compared with the
experimental values.  I could find no indication as to whether the
reported calculations used the UHF or ROHF approaches.  This is
important because the calculated heats for each approach produces
significantly different results.  I consequently repeated the
calculations on the 6 hydrocarbon radicals given in the above paper
and found the following heats of formation (kcal mol^{-1}).  I also
give the values of S^2 and the reported heats as well as the current
experimental values (from J. Phys. Chem. Ref. Data, 17 (1988) Sup. 1)
for comparison below.

           AM1           AM1    Reported
           UHF    S^2    ROHF      AM1     Exp.     Exp. - Calc.(ROHF)
CH3       29.95  0.7613  31.25    31.25   34.8(3)       3.6
C2H3      60.46  0.8589  64.78    64.78   63.4(10)     -1.4
C2H5      15.49  0.7619  18.14    15.48   28           10
CH2CHCH2  30.20  0.9300  38.58    38.58   39            0
(CH3)2CH   3.61  0.7622   6.80    10.07   22.3(6)      15.5
(CH3)3C   -6.14  0.7623  -2.66    -2.66   11.0(6)      13.7

For C2H5 the authors seem to have reported the UHF value.  I was not
able to reproduce the reported result for (CH3)2CH with either the
ROHF or the UHF result.  I was able to reproduce the calculated heat
of formation for this species given by Dewar, M.J.S. and Thiel, W., J.
Am. Chem. Soc., 99 (1977) 4907, using the MNDO Hamiltonian.  In this
earlier work the authors explicitly stated that they used ROHF heats
of formation calculated with MNDO.  It is noted that for the remaining
species the reported AM1 heats of formation also appear to be the ROHF
values.  While agreement with the experimental values are not terrific
and appear particularly bad for the larger saturated hydrocarbon
species, they are adequate for my purposes given that these errors can
be taken as typical for calculated heats of formation of analogous but
larger open shell hydrocarbon radicals.

Does anybody know of any further studies, using the AM1 Hamiltonian,
on radicals where the calculated heats of formation have been compared
with experiment?  It seems to me, from the above table, that the UHF
results are not useful in comparison with the ROHF results, and that
the extent of spin contamination tends to increase with molecular size
and unsaturation, being close to 0.75 for the most saturated species.
For instance, I repeated the calculation on triplet linear C19 given
by Novoa et al., Inorganica Chemica Acta 198-200 (1992) 133, using UHF
and ROHF.  I obtained the same heat of formation as these authors with
no assumed symmetry in the linear chain and using the UHF AM1
Hamiltonian.  However, I obtained an S^2 of 4.356 (should be 2) and a
ROHF heat of formation of 678.72 kcal mol^{-1}, which is greater than
the UHF value by 24.17 kcal mol^{-1}!  I would very much appreciate
any helpful comments or literature references which deals with any of
the above issues I've raised here about AM1 and its applicability of
open shell systems.

Regards,

Ryan Bettens,
OSU Physics Department,
[email-masked]

                       ________________________
----------------------| Responses to Questions |----------------------
                       ------------------------

From:  Andrew Holder,
       University of Missouri,              
       Department of Chemistry               

   You have indeed found some the reporting errors in the original AM1
paper.  The problem with your larger systems is that the description
of an open shell system of this size will usually require more than 
one determinant.  (That is what the bad spin contamination values 
indicate.)  A better way to handle this is to do everything at the 
same level of configuration interaction (CI).  Be sure that you 
compare everything in a reaction profile or an analogous series 
computed at the same level of theory.

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

From:  John M. McKelvey,
       Eastman Kodak Company, Rochester NY              
       Computational Science Laboratory               

There is no proper ROHF in QCPE MOPAC or AMPAC codes!  It is a half-
electron method supplemented by a 3x3 CI.  I THINK there is proper 
ROHF for AM1 and PM3 and MNDO in GAMESS.

             ----------------------------------------
             
Comment to J.M. McKelvey's response.

It seems I had too loosely used the ROHF terminology in my posting.  
When I stated that Dewar and Thiel, ". . . explicitly stated that they
used ROHF heats of formation calculated with MNDO" I ignorantly 
mistook the statement actually made by the authors about the method 
they used, i.e., the half-electron method, for ROHF (which from what I
gather from various modern books on ab initio theory is Roothaan's 
open shell method).  Therefore, in the above posting wherever ROHF 
appears, read: half-electron method.

John McKelvey pointed me toward the publication:  Dewar, M.J.S. and 
S. Olivella, J. Chem. Soc. Faraday Trans. II 75 (1979) 829, which 
compared the calculated heats of formation and geometries for the 
half-electron (HE) method with what is now meant by ROHF using 
MINDO/3.  Their results showed that the ROHF heats of formation were 
"systematically more negative than the HE ones".  They also found that
for the molecules studied the r.m.s. differences in the heats of 
formation for open shell doublets, triplets and singlets were 2.2, 2.5
and 3.7 kcal mol^{-1} respectively.

For those interested, in the ROHF method the wavefunction is (a) 
variationally optimized and (b) is an eigenfunction of the S^2 
operator, unlike UHF wavefunctions.  In the HE method the wavefunction
is not variationally optimized and is the same for analogous open 
shell states of different multiplicity.  A non-mathematical 
description of the HE method can be found in the introduction of the 
Dewar and Olivella paper, above.

************************************

Date: Tue, 08 Nov 1994 15:56:15 -0400 (EDT)
Subject: Summary: UHF vs ROHF in semiempirical calc's of radical cations
To: [email-masked]
Message-id: <[email-masked]>
X-VMS-To: IN%"[email-masked]"
X-VMS-Cc: ZUILHOF
MIME-version: 1.0
Content-type: TEXT/PLAIN; CHARSET=US-ASCII
Content-transfer-encoding: 7BIT

Dear CCl'ers,
As promised I hereby send you the summary of all responses I received on my 
questions about the use of ROHF versus UHF in semiempirical calculations of
radical cations. This summary starts with the original posting, followed by the
responses.
A big THANKS to all those who responded!

                                              Han Zuilhof

##############################################################################
Original posting:

Subj:	ROHF versus UHF in semiempirical calc.'s of radical cations

Dear CCl'ers,

Semiempirical programs such as MOPAC allow the properties of radicals
and radical ions to be calculated with both the Restricted Open-shell
Hartree-Fock method (ROHF) and the Unrestricted Hartree-Fock (UHF) methods.
I want to calculate reactionpaths for nucleophilic attack on a series of
radical cations.

Does anyone know of 
1) any apriori reasons why one method might be preferred over the other in 
such calculations?
2) literature data in which the performance of the ROHF and UHF methods 
for the study of radical cations is compared directly?

Please report directly to me, and I'll summarize to the net.
Thanks in advance,
                    Han Zuilhof

##############################################################################

A very informative answer, dealing with many issues came from 
prof. Thomas Bally:

From:	IN%"[email-masked]"  2-NOV-1994 08:02:49.67
To:	IN%"[email-masked]"
CC:	
Subj:	RHF vs. UHF

Dear Mr. Zuilhof,

You were asking on CCL about using UHF vs.RHF for calculating the attack of
nucleophiles on radical cations.

The question of UHF vs. RHF is a thorny one and it becomes thornier when
you go to semiempirical methods, even if you disregard for the moment the
problems posed by spin contamination which I will address below.

By allowing alpha and beta-electrons to occupy different spatial MO's (which is
the essence of the UHF model) you take into account a special form of electron
correlation called "spin polarization" which expresses itself experimentally in
negative hyperfine coupling constants in ESR spectra and is hence a physical
phenomenon. However, the UHF model has a tendency to overestimate the extent of
spin polarization which results in an overestimation of the stability of
systems where spin polarization is important (for example systems containing
allylic radical moieties such as they occur quite often in radical ions) if
compared to systems where this effect is less important.

In semiempirical methods, electron correlation effects are already "absorbed"
somehow in the parameters and therefore, in systems with high spin polariza-
tion, you count some of these twice in semiempirical UHF calculations. As a
consequence, the stability of systems where spin polarization is important (see
above) is often *absurdly* overestimated by semiempirical methods.

On the other hand, if you do RHF (I suppose you are referring to the "half-
electron" method which is used rather than the "real" Roothaan ROHF method in
Dewar's modelsbut this makes no difference), you artificially suppress spin
polarization entirely (it cannot occur if alpha and bete-electrons are
constrained to occupy pairwise the same spatial MO's). Hence the stability of
systems where this is important (see above) is *underestimated* with this
model.

To this you have to add the problem of "high-spin contamination" in UHF
wavefunctions which are not eigenfunctions of the S**2 operator and hence
cannot be classified as singlets, doublets, triplets. Most decent programs will
tell you about the *expectation value* of a UHF wavefunction with regard to
S**2 (usually called <S**2>) and the deviation of this expectation value from
the value for a pure singlet(=0), doublet (=0.75), triplet (=2.75) etc. tells
you how far you have gone astray from pure spin states in your calculation. As
a rule of thumb, a UHF wavefunction is acceptable if <S**2> deviates less than
0.1 from the "correct" value. You will find that this deviation is often
exceeded, especially in systems with high spin polarization.

So what can you do? Unfortunately I cannot give you a good answer. Due to the
above-mentioned problems and some other shortcomings of semiempirical methods
(problems with small rings) I have personally all but given up semiempirical
calculations on radical ions and switched to ab-initio based methods which have
become quite affordable with the advent of modern workstations. On the other
hand you may get away with semiempirical methods for large systems if

  (a)  You do have cases where spin polarization is unimportant
       (i.e. if the difference between UHF and RHF heats of formation
       is, say, less than 5 kcal/mol)

  (b)  If the <S**2> values for doublets are <0.85.

The grist of the matter is, that you should always do *both* RHF and UHF
calculations. If the results do not differ substantially, both with regard to
geometries and relative energies, you are probably safe (unless you have small
rings, but that's another can of worms).

Don't hesitate to get back to me if you have follow-up questions or if you need
guidance for getting started with ab-initio calculations (of such calculations
are feasible for your systems)..

  thomas bally

P.S. please don't reply to the address from where you got this message (I use
     it exlusively for CCL) but to the one indicated below:

------------------------------------------------------------------------------
|  Prof. Thomas Bally                    |  E-mail:  [email-masked]   |
|  Institute for Physical Chemistry      |                                   |
|  University of Fribourg                |  Tel:     011-41-37 826 489       |
|  Perolles                              |  FAX:     011-41-37 826 488       |
|  CH-1700 FRIBOURG                      |                                   |
|  Switzerland                           |                                   |
------------------------------------------------------------------------------

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

Both Serge Pachkovsky and John McKelvey commented that MOPAC/AMPAC etc do not
implement the REAL/CORRECT ROHF method, but rather Dewar's half-electron method:

From:	IN%"[email-masked]"  2-NOV-1994 03:58:29.52
To:	IN%"[email-masked]"
CC:	
Subj:	RE: CCL:ROHF versus UHF in semiempirical calc.'s of radical cations

Dear Han,

In your message to the CCL you write:
> Semiempirical programs such as MOPAC allow the properties of radicals
> and radical ions to be calculated with both the Restricted Open-shell
> Hartree-Fock method (ROHF) and the Unrestricted Hartree-Fock (UHF) methods.
                       ^^^^

This is not exactly true. Semiempirical programs usually implement half-electron
method (and Mopac is not an exception), not ROHF, which is notorious for the
poor SCF convergence properties. Although half-electron method often gives
results one would expect from the exact ROHF treatment (and in some cases is
exactly equivalent to it), there is one important distinction: electronic
energy computed within the half-electron approximation is not variational.

With my best regards,

Serge Pachkovsky.

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

Careful...NONE of the semi-empirical methods derived from MOPAC or AMPAC does 
the CORRECT ROHF. They all use a half-electron method of Dewar followed by a 
small C.I. for clean up. Gaussian and GAMESS will do the ROHF for the 
semiempirical methods correctly, I believe.

-- 

John M. McKelvey			email: [email-masked]
Computational Science Laboratory	phone: (716) 477-3335
2nd Floor, Bldg 83, RL
Eastman Kodak Company			
Rochester, NY 14650-2216

---------------------------------------------------------------------------
the matter of spin contamination occurring in UHF calculations was adressed
by several people:
Dave Ewing writes:

From:	IN%"[email-masked]"  "DAVID W. EWING (216) 397-4742"
To:	IN%"[email-masked]"
CC:	
Subj:	RE: CCL:ROHF versus UHF in semiempirical calc.'s of radical cations

With ROHF you don't have to worry about spin contamination.

Dave Ewing
John Carroll University
[email-masked]
----------------------------------------------------------------------------

and Irene Newhouse wrote:

From:	IN%"[email-masked]"  3-NOV-1994 12:35:15.41
To:	IN%"[email-masked]"
CC:	
Subj:	RE: CCL:ROHF versus UHF in semiempirical calc.'s of radical cations

There are people who are VERY suspicious of ROHF -- people like Andy Holder,
one of the authors of AMPAC v4.  On the other hand, UHF can, especially very
close to the transition state, suffer from 'spin contamination', that is,
states with the wrong spin mix in.  At very stage, SZ**2 is plotted, so that
you can TELL if it's happening.  Which way you do it, depends on if you're
trying to qualitatively describe experiments, or do a thorough theoretical`
investigation.  If the former, I'd try it both ways & be happy if the results
are about the same.  If the latter, all either method is good for is a guess
to plug into ab initio computations!

Good luck!`
Irene Newhouse 
###############################################################################
###############################################################################

******************************************************************************
**   Dr. Han Zuilhof           **  e-mail: [email-masked]  **
**   Department of Chemistry   **  (optional: [email-masked]) **
**   University of Rochester   **                                           **
**   Rochester, NY, 14627      **  Fax:   (716) 473-6889                    **
**   USA                       **  Voice: (716) 275-2219                    **
******************************************************************************
**                                                                          **
**                        "Excite a photochemist!"                          **
**                                                                          **
******************************************************************************

>From [email-masked]  Mon Nov 14 09:37:32 1994
From: <[email-masked]>
Date: Mon, 14 Nov 1994 09:11:33 -0500 (EST)
Subject: Spin contamination, effect on energy and structure.
To: [email-masked]

Dear Netters,

I have a number of questions regarding the effects of spin=20
contamination on total electronic energies and structures.  My=20
understanding of spin contamination is that unrestricted Hartree-Fock=
=20
(UHF) wave functions are not eigenfunctions of the total spin=20
operator, so the electronic wave function of interest can be=20
contaminated by functions corresponding to states of higher spin=20
multiplicity.  This brings me to my questions:

1.  Given that, (a) we have performed an ab initio study on an open=
=20
shell molecule using a single-determinant wave function (i.e.,=20
variational), (b) we have found its theoretical equilibrium geometry=
=20
for the ground electronic state, and (c) have projected out ALL=20
contaminating higher spin multiplicity states along the way to the=
=20
minimum energy geometric configuration.  Will the total calculated=
=20
electronic energy be the lowest possible calculated energy for the=
=20
given basis set and Hamiltonian?

2.  What can be said about 1, above, regarding the total calculated=
=20
electronic energy when a perturbation treatment to the configuration=
=20
interaction is introduced, e.g., MP4(SDTQ)?

3.  What can be said regarding the total calculated electronic energy=
,=20
if in the case of 2, above, condition (1c) is not fully met, i.e.,=
=20
some spin contamination remains while going to the optimum geometry.

4.  Given that different spin states correspond to different=20
equilibrium geometries, is an optimized structure, which had=20
significant spin contamination in its electronic wave function all th=
e=20
way down to its =D2minimum=D3 energy geometrical configuration, some =
kind=20
of mixture of structures involving the state of interest and the=20
different contaminating higher spin multiplicity states?

5.  Regarding semiempirical calculations, in the paper of Novoa et=
=20
al., Inorganica Chemica Acta, Vol. 198-200 (1992) 133, the heats of=
=20
formation of some very large carbon clusters were calculated.  In=
=20
their paper the authors state:  =D2For an odd-membered linear C_n wit=
h=20
13 <=3D n < 20, the AM1 calculations predict the triplet state to be=
=20
more stable than the singlet, due mainly to the spin contamination of=
=20
the UHF calculations.=D3  The calculated stabilities of these larger=
=20
linear carbon clusters are not expected based on what is known for th=
e=20
smaller linear clusters where, for odd n (n > 1), the singlet states=
=20
are more stable than the triplet states.  My question is this, is it=
=20
possible that the authors are not correct regarding their reason for=
=20
the for greater stabilities of the triplet states?  (It is not my=
=20
intention to attack the above authors work.  I merely wish to evaluat=
e=20
the quoted heats of formation because I require some kind of estimate=
=20
for the heats of formation for these species.)

Regards,

Ryan Bettens,
OSU Physics Department,
[email-masked]