Spin contamination & AM1 "ROHF" versus UHF



 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,
 BETTENS (+ at +) MPS.OHIO-STATE.edu
                        ________________________
 ----------------------| 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,
 BETTENS (+ at +) MPS.OHIO-STATE.edu
                        ________________________
 ----------------------| 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.