Subject: Summary for "Calculation on biradical"

From: <[email-masked]>
Date: Tue, 29 Aug 1995 15:33:48 +0200 (MST)

Dear netters,

 some days ago I posted a question about calculations on a biradical. The 
original posting and the answers are attached. 

 Thanks to all who responded!
 
 Thorsten Koch

The original question:

Dear netters,

 At the moment I am doing ab initio calculations on 
1,5-Dehydronaphthalin, so it is supposed to be a biradical. It is an 
isomer to a ordinary closed-shell molecule which I am interested in. So 
to compare the energies of the molecules, I am interested in the 
singlet-groundstate of the biradical.
 But this doesn't seem to be a trivial problem. It was no problem at all 
to get the wavefunction for the triplet state. So I optimized the 
geometry for the triplet state and tried to get some singlet 
wavefunctions for this geometry.

 Here is what I got:

 1. UHF, 6-31G*, mult=3 : E=-382.059168, <S**2>=2.2470 (2.0346 after 
annihilation of first spin contaminant)

 2. UHF, 6-31G*, mult=1 : E=-381.906022, <S**2>=0.0000
 This result was the same as the result with RHF. Especially the spin 
desity was zero everywhere. So I thought that this wouldn't be a proper 
description of a biradical.

 3. UHF, 6-31G*, mult=1, GUESS=mix
 This calculation didn't converge at all.
 
 4. In fact I was able to get another singlet state wavefunction. The 
only way to get it was from a stability analysis and optimization from the 
wavefunction obtained with 2. The result was: E=-382.086893, 
<S**2>=2.2807 (8.1431 after annihilation of the first spin contaminant)
At least the spin density gave two unpaired electrons with different spin 
at the right carbons.

 Because dft is known to give lower spin contaminations, I tried it this way:

 5. UHF, 6-31G*/B3LYP, mult=3 : E=-384.510256, <S**2>=2.0136 (1.9996)
 
 6. UHF, 6-31G*/B3LYP, mult=1 : E=-384.484437, <S**2>=0.000
 Again the same result as with RHF, spin density zero everywhere.

 7. UHF, 6-31G*/B3LYP, mult=1, GUESS=mix : E=-384.518499, <S**2>=0.9392 
(0.1510 after annihilation) This calculation had no difficulties to converge.

 ROHF calculations for the singlet states always gave the same 
wavefunction as for RHF calculations.

 So there are a couple of questions:
 
 a) Why is it so difficult to get a UHF singlet wavefunktion for a normal 
biradical? 

 b) Can I rely on the energies obtained when the spin contamination is 
about 1.0 (Point 7.)? This is especially problematic because I want to 
compare it to the energiy of a closed-shell isomer. Since I got this 
results with 6-31G*/B3LYP it seems to be the best I can get...
Is a geometry optimization sensible in this case?

 c) What can I say about the wavefunction of point 4.? The spin 
contamination is about 2.3 before annihilation (not a very small value 
for a singlet ;-) and about 8.1 AFTER annihilation of the first spin 
contaminant. I have no idea what this should mean!

 d) Am I right when I think that the spin density of a system with to 
unpaired electrons in the sigma frame separated far enough from 
another should exhibit these two electrons at the carbons they are 
attached to?

e) Having all these problems in mind: Which way should I go to determine 
the energiy of the singlet ground state?

Sorry for this somewhat lengthy question, but I'll summarize the 
responses, when there is interest :-)
-----------------------------------------------------------------------

The answers:

From: Steve Gwaltney <[email-masked]>

Your molecule sounds like an open-shell singlet.  If it is, you need
two determinants to describe the wavefunction.  Thus, any Hartree-Fock
calculation will fail.  What you need is a MCSCF calculation, or some
other way to deal with two determinantal wavefunction.

Steve

From: Christopher J Cramer <[email-masked]>

>  a) Why is it so difficult to get a UHF singlet wavefunktion for a normal 
> biradical? 
> 
   Because an open-shell singlet is a two-configuration wavefunction. It
requires an MCSCF approach to do it properly. The S^2=1 calculations are for
so-called 50:50 wavefunctions (half singlet, half triplet).

>  b) Can I rely on the energies obtained when the spin contamination is 
> about 1.0 (Point 7.)? This is especially problematic because I want to 
> compare it to the energiy of a closed-shell isomer. Since I got this 
> results with 6-31G*/B3LYP it seems to be the best I can get...
> Is a geometry optimization sensible in this case?
> 
   No, and no. You can adopt the so-called sum method of Ziegler and
double-the energy difference between the triplet and the 50:50 wavefunction
to estimate the energy of the open-shell singlet. Although a terrible
approach at the HF level, comparing UDFT to RDFT seems to be OK for S-T gaps
(we have several papers on this). However, your B3LYP functional includes HF
exchange, which will degrade this comparison.

>  c) What can I say about the wavefunction of point 4.? The spin 
> contamination is about 2.3 before annihilation (not a very small value 
> for a singlet ;-) and about 8.1 AFTER annihilation of the first spin 
> contaminant. I have no idea what this should mean!
> 
   It's garbage.

>  d) Am I right when I think that the spin density of a system with to 
> unpaired electrons in the sigma frame separated far enough from 
> another should exhibit these two electrons at the carbons they are 
> attached to?
> 
   For the closed-shell singlet, there is probably hybridization and orbital
splitting. For the open-shell and the triplet, your intuition is right. Much
of the literature on 1,4-didehydrobenzene discusses these topics.

> e) Having all these problems in mind: Which way should I go to determine 
> the energiy of the singlet ground state?
> 
   MCSCF is the best approach. Your DFT results are quite curious because
they indicate the open-shell singlet to be lower in energy than the triplet.
This is pretty unusual, and I'm not sure I believe it.

   If you would be interested in seeing some of our work on these issues, I'd
be happy to send you some papers.

Best regards,

Chris

-- 

From: Christian Koelle <[email-masked]>

> 
> e) Having all these problems in mind: Which way should I go to determine 
> the energiy of the singlet ground state?
> 
Zwei Wege: einer gut der andere besser:

1) CI Rechnung. Bei einem Diradikal sollte als aktiver Raum HOMO und LUMO
ausreichen. Wir benutzen bei unseren Rechnungen mit einer semiempirischen 
Methode die HOMO-LUMO Einfach- und Doppelanregung, wobei die 
Doppelanregung fuer den Singulett entscheidend ist.

2) MCSCF-Rechnung. Auch hier HOMO und LUMO als aktiver Raum.

Hoffe geholfen zu haben.

Gruss

Christian Koelle

From: Hans Ueli Suter <[email-masked]>

Lieber Herr Koch,

ich wuerde Ihr Problem unter einem etwas anderen Gesichtspunkt sehen:
Was Sie effektiv rechnen sind die angeregten Zustaende von ...
Bei einem Programm wie Gaussian, das keine direkte Symmetriebindung der
Orbitale kennt ist es natuerlich schwierig diese in der noetigen
Weise zu besetzen. Wie dem auch sei, das Keyword wuerde lauten
Guess=Alter und dann am ende des Inputs waeren die Numern vom HOMO und
vom LUMO einzusetzen, dann wuerde er eben das HOMO einfach besetzen und
das LUMO einfach, womit Sie das hoffentlich richtige Biradikal haetten.
Ich war nie in der Lage das auch fuer Dichtefunktionale durchzufuehren,
aber UHF sollte gehen. Andererseits koennen Sie versuchen mit CIS
die angeregten Zustaende zu rechnen (Das Biradikal wird auch darunter
sein). Vom Niveau her ist "CIS" etwa vergleichbar mit HF und nicht mit CI.
Die Resultate sind demzufolge mit Vorsicht zu geniessen. Hier wuerde jetzt
natuerlich der Hinweis folgen, dass man angeregte Zustaende besser mit ..., 
aber ich spare mir das, sollte eh klar sein. Ich hoffe, dass Ihnen das hilft.

                                   mit freundlichen Gruessen
                                        H.U. Suter

From: Frank Jensen <[email-masked]>

	Thorsten,
	UHF wave functions that are spin contaminated to the
degree you describe are essentially useless, even your triplet
UHF is pretty bad. The energetics will be all garbage. You need
a small MCSCF (or equivalently GVB) to get a reasonable zero'th
order description. For your singlet biradical a [2,2]-CASSCF
may be OK (or an OSS GVB type wave function which is the same).
But you won't get any dynamic correlation. To include that you
need a MRCI, or perhaps a MP2 on top of the MCSCF.
	Frank

From: [email-masked]

Singlet biradicals are no trivial problem. You see, your calculations
are on the *lowest* singlet, not the biradical. The calculations are 
converging to (what I believe is called) a Renner-Teller distorted 
configuration of the orbitals because it wants to do this

                                                       ------
        U         D
     -------    ------       distorts to
                                                         UD
                                                       ------

What you want to do is not UHF, where your spins are all crazy mixed up. What 
you want to do is CI, where you can specify exactly what sort of configuration 
you want, and send it on it's merry way. I believe QCISD and CASSCF are the 
methods du jour for these sort of problems.

matt