Subject: SUMMARY: g98: TD-DFT geometry optimizations ?

From: Johannes Weber <[email-masked]>
Date: Thu, 11 Apr 2002 13:54:35 +0200 (MET DST)

Dear CCL members,

Two weeks ago I posted a g98 related question concerning geometry
optimizations and frequency calculations with the time-dependent DFT
method. Here is the original question:
-------------------------------------------------------------
Dear CCL members,

I have a Gaussian 98 related question:

According to the g98 users manual (p. 155, TD keyword) it should be
possible to do Time Dependent DFT geometry optimizations of the N-th
excited state using the Root=N option. 

However, I can´t find a proper input. The easiest input line

 #P TD=(Root=1) B3LYP/STO-3G Test DENSITY OPT=Z-Matrix

leads to an error in link 1002:

 (Enter /vol/chemie/libexec/g98.A9/g98/l1002.exe)
 Minotr:  Closed-shell wavefunction.
 Cannot mix post-SCF and DFT.
 Error termination via Lnk1e in /vol/chemie/libexec/g98.A9/g98/l1002.exe.
 
I played around with several other options, tried to do the job in a two
step procedure, like it is described for the CIS method in
Forsman/Frisch´s "Exploring Chemistry with Electronic Structure Methods",
all without success. Is it possible to do these calculations?  Anybody out
there with experience in this topic? The ultimate goal would be to perform
a normal mode frequency calculation for excited states within the TD-DFT
formalism. 

Regards, :-)ohannes.
---------------------------------------------------------------------

I got several answers, my special thanks go to Dr. Doug Fox and Prof.
Galina Orlova for the fruitful discussion. Other people which I´m indebted
to are Wolfgang Roth and Rainer Remenyi. To put it in a nutshell: 

Yes, you can do TDDFT geometry optimizations and frequency calculations. 
However, since no gradients and higher derivatives are available for this
method in g98, you have to resort to numerical optimizations, namely using
the EF-algorithm with the opt=(Enonly,Z-Matrix) and freq=(Enonly) 
keyword. The optimization has to be performed in Z-Matrix coordinates
and you need an initial guess for the force constants, which Dr. Fox
recommended to take from a CIS calculation (optimization and frequency
calculation). You can also try an take the force constants out of a
different calculation with the RCFC option in the opt keyword. Opt and
freq calculation have to be done in subsequent steps, a one-step procedure
with opt and freq in the same input card will fail.

To my (very limited) experience the usability of TDDFT optimizations is
restricted by the fact that the numerical derivatives suffer from long run
times (especially the freq calculations are really expensive for
polyatomic molecules!) and a quite low accuracy (this is at least what I
found from in the calculated frequencies). Until now, I wasn´t able to
use force constants from other than CIS calculations for the initial
guess.

Best regards and thanks again to all people who answered. The relevant
original answers are appended below. 
  
 :-)ohannes. 

======================================================================
answer 1:
From: Doug Fox <[email-masked]>
Subject: Re: CCL:g98: TD-DFT geometry optimizations ?

  Dr Weber,

   The manual is in error and there is not an implementation of
analytic gradients for TD-DFT nor is there for TD-HF (aka RPA).

   You can get various roots and if you use OPT=(RCFC,EnOnly), taking the
force constants perhaps from a single point CIS FREQ calculation, you
could do a geometry optimization.  Note that this will be more practical
for smaller systems. 
======================================================================

answer 2: Prof Orlova forwarded an email to me that she had received
earlier from Jim Hess, [email-masked]. 

From: Galina Orlova <[email-masked]>
Subject: Fwd: Re: TDDFT (from WWW)

----- Forwarded message from [email-masked] -----
    Date: Mon, 25 Mar 2002 11:49:06 -0500
    From: [email-masked]
Reply-To: [email-masked]
 Subject: Re: TDDFT (from WWW)
      To: [email-masked]

Dear Galina,
      Since there aren't analytic gradients available for TD methods, 
normal geometry optimizations are not possible with this keyword. The 
error you get with TD-DFT occurs when Gaussian actually tries to 
calculate the analytic gradient, but fails to find the necessary data. 
The only way to do a td-dft optimization is with numerical gradients. To

do this, use "opt=(enonly,rcfc)". This option requires that the 
optimization be done in symbolic z-matrix coordinates and that initial 
force constants be read from the checkpoint file or from cards. The 
force constants would typically come from an analytic frequency. In this

case, CIS frequencies seem to be the logical choice. As for your TD-HF 
optimizations, I'm afraid to say that they are not truly TD-HF 
optimizations. The code should give an error message, but instead it 
does the optimization using something very close to a CIS gradient. 
Since the optimizations you already completed were done with gradients 
based upon CIS, if you choose to go ahead with the numerical TD-HF 
geometry optimizations, these final geometries should be fairly close to

the TD-HF minimum so I'd recommend using them as initial structures.

Regards,
Jim Hess
---------------------------------------------------------------------
...
Dr. Galina Orlova
Vis. Res. Professor
Department of Chemistry and Biochemistry
University of Guelph
Guelph, ON, Canada N1G 2W1
tel.  work  (Int+1)(519)824-4120 ext 3102
      home  (Int+1)(519)829-3119
fax                (519)766-1499 
e-mail  [email-masked]
=====================================================================

answer 3:
From: Wolfgang Roth <[email-masked]>
Subject: Re: CCL:g98: TD-DFT geometry optimizations ? 

Hallo Johannes,

ich habe zwar schon eine ganze Weile nichts mehr mit Gaussian gerechnet
und kenne auch die neuesten Features nicht, ABER:

Zu dem von Dir beschriebenen Problem
>According to the g98 users manual (p. 155, TD keyword) it should be
>possible to do Time Dependent DFT geometry optimizations of the N-th
>excited state using the Root=N option. 
kann ich Dir folgendes sagen:
Der TD-DFT output von G98 Rev a.7 und a.9 sieht genauso aus wie die
AUsgabe, die man bei einer CIS-Rechnung bekommt. Deswegen gehe ich mal
davon aus, daß die Programmierer die Algorithmen wiederverwendet haben.
Das "Root" Keyword beschreibt nur den "state of interest" d.h. den
gewünchten elektronischen Zustand, für den die gewählte Rechnung
durchzuführen ist. Ob man oder ob man nicht Geometrieoptimierungen
durchführen kann ist bei der Berechnung der Wellenfunktion nicht von
belang. Meines Wissens gibt es in keinem Programm einen Algorithmus zur
Geometrieoptimierung mit Time dependent Methoden, auch nicht für TD-HF,
das
ja eigentlich der CIS-Methode sehr ähnlich/vergleichbar sein sollte.
Das Problem ist, daß es wohl keine Gradienten für den TD (EOM??)
Algorithmus gibt.

> Is it possible to do these calculations?  
Ich denke: NEIN
ich habs auch mal über die numerische Gradienten versucht, da tut Gaussian
dann ywar ne ganze Weile (3N TD-DFT single points lang) was, steigt dann
aber mit Fehlermeldung/ohne Ergebnis aus.

>Anybody out there with experience in this topic? The ultimate goal would
be to perform
>a normal mode frequency calculation for excited states within the TD-DFT
>formalism. 
Ja, das wär wirklich nett, v.a. weil CASSCF Freqs mit Gaussian mühsam
sind.. 

Schöne Grüße
Wolfgang

==W=R===========================W=o=l=f=g=a=n=g==R=o=t=h==
  TU Brauschweig
  Institut f. Physikalische und Theoretische Chemie
  Hans-Sommer-Strasse 10        Tel ++49 (0)531 391 5346
  38106 Braunschweig            Fax ++49 (0)531 391 5396
======== http://www-public.tu-bs.de:8080/~wroth/ =========

=====================================================================

answer 4:

From: Rainer Remenyi <[email-masked]>
Subject: Re: CCL:g98: TD-DFT geometry optimizations ? 

Hi Johannes,

Try  UB3LYP instead of B3LYP.

Regards,

Rainer Remenyi
Universität Heidelberg
=====================================================================

answer 5: The second answer of Prof. Orlova, after I asked her if she also  
          calculated the frequencies.

Hi, Johannes

I guess, if G98 does not have analitic gradients for TD then it does not
have analitic second deriv. either. So you must run frequencies
numerically (freq=numer). If you did do so and it still doesn't work you
better ask the gaussian helpers. Unfortunately, my systems are so large
that even numerical optimization is a great problem. Thus I don't have
experience with numerical frequencies within TDDFT. 

Jim means the TD-hartree-Fock (not TDDFT!) optimization I've run
"successfully" analitically. In fact, G98 must abort this job. For some
reason it was running but Jim's guess is that in fact G98 performed the
CIS-HF optimization instead of TD-HF. Actually,the CIS-HF results might be
very close to those from TD-HF optimization with numerical derivatives. 

best regards, Galina

Dr. Galina Orlova
Vis. Res. Professor
Department of Chemistry and Biochemistry
University of Guelph
Guelph, ON, Canada N1G 2W1
tel.  work  (Int+1)(519)824-4120 ext 3102
      home  (Int+1)(519)829-3119
fax                (519)766-1499 
e-mail  [email-masked]
========================================================================

answer 6: The second answer of Dr. Fox, after I tried out a numerical 
          frequency optimization -- with peculiar results. 

From: "Cust. Service Doug" <[email-masked]>
Subject: Re: CCL:g98: TD-DFT geometry optimizations ?

  Dr. Weber,

   There is no option to do numerical or enonly frequencies as a
part of an optimization.  
...
  As to your first point clearly the attempt a multi-step job, ie. OPT
FREQ, is generating a bad keyword.  As this is still present in the
current revision I can only suggest you do this as a two step job,
OPT=(RCFC,Enonly) as the first and then FREQ=Enonly for the second.
I will have the combination examined for a future release but we don't
plan any further bug fixes for G98 at this point.  Note that it has
created the bad combination RTD-B3LYP.  It may also have a problem 
with correctly generating FREQ=EnOnly but I can't confirm that because
it fails before actually generating that route.

   Numerical frequencies do suffer from much lower precision.  The
step size is a compromise between the infinitesmal step which 
would be the closes approximation to the analytic derivative but
which has such a small energy difference that it is numerically
inaccurate or a larger step which gives a numerically stable
energy difference but might step outside of the harmonic region
and introduce errors that way.  In the case of very low frequencies
you are probably better off using a somewhat larger step as the 
low frequencies will require a larger step to get a significant
energy difference.  Smaller steps might help the higher frequencies
but that is less important for the optimization.
...

> 
> Dear Dr. Fox,=20
> 
> thanks again for your quick response.  After the easter vacation I now
> come back to the TDDFT geometry optimization and frequency calculation
> problem and pester you with further questions. Please, be patient with
me!=
> =20
> :-)  I tested several cases with a minimal basis on para-benzoquinone.
My
> results are as follows:=20
> 
> (1) I managed to do a numerical frequency calculation after the
> optimization job. However, it seems not possible to do a geometry
> optimization and a frequency calculation in a single step! Having done a
> CIS-optimization and freq calculation in advance the input card
> 
> > # TD=3D(Root=3D1,Read) RB3LYP/STO-3G INTEGRAL(GRID=3DULTRAFINE) Test
> > DENSITY=3Dcurrent Geom=3DCheck Guess=3DRead OPt=3D(EnOnly,RCFC)
Freq=3DEN=
> Only  =20
> 
> leads - after a sucessful TD-DFT geometry optimization - to the
following
> error message in the freq run:
> 
> >...
> > ------------------------------------------------------------
> > #P Geom=3DAllCheck Guess=3DTCheck Test RTD-B3LYP(FC)/STO-3G Freq
> > ------------------------------------------------------------
> >  QPERR --- A SYNTAX ERROR WAS DETECTED IN THE INPUT LINE.
> >  GUESS=3DTCHECK TEST RTD-B3LYP(FC)/STO-3G=20
> >                             '
> > Last state=3D"GCL"
> > TCursr=3D  820 LCursr=3D  44
> > Error termination via Lnk1e in /vol/chemie/libexec/g98.A9/g98/l1.exe.
> > Job cpu time:  0 days  0 hours  0 minutes  0.1 seconds.
> > File lengths (MBytes):  RWF=3D  105 Int=3D    0 D2E=3D    0 Chk=3D
12 S=
> cr=3D1
> > Segmentation Fault
> >...
> 
> In contrast,  the following two step procedure seems to work:
> 
> > %chk=3Dpbq.chk
> > # TD=3D(Root=3D1,Read) RB3LYP/STO-3G INTEGRAL(GRID=3DULTRAFINE) Test
> >    DENSITY=3Dall Geom=3DAllcheck Guess=3DRead OPt=3D(EnOnly,RCFC)=20
> >
> >
> Why is this the case? Does the single step procedure always assume
> analytical frequencies, even though Freq=3DENOnly was specified
originally?=
> =20
...
======================================================================

----------------------------------------------------------------------
| Johannes Weber                     | Email:                        |
| Universitaet zu Koeln              | [email-masked]   |
| Institut fuer Physikalische Chemie |                               |  
| Lehrstuhl Prof. Dr. G. Hohlneicher | tel: 0049-(0)221-470-4812     |
| Luxemburger Str. 116               | fax: 0049-(0)221-470-5144     |
| 50939 Koeln                        |                               |
----------------------------------------------------------------------