Subject: SUMMARY: DFT B3LYP/6-31+G(d,p) Convergence Problem

From: "M. Nicklaus" <[email-masked]>
Date: Fri, 26 Jun 1998 20:19:48 -0400 (EDT)

Dear CCL'ers,

This is the summary of the responses (slightly edited for brevity)
I received to my question.  The stepwise approach of going to
successively higher basis sets, suggested by several people,
worked just fine.  Thanks to all who responded.

Marc

================= Original Question (Tue, 9 Jun 1998) =================

After optimizing a structure (60 atoms, 37 non-H) at the HF/3-21G* level,
I'm trying to calculate a DFT single point energy at the B3LYP/6-31+G(d,p)
level, reading geometry and initial wave function from the checkpoint file.

After nearly seven days (!) of CPU time on a very fast machine, the job
aborts without having converged, with the following end of the output file:

[snip]
 >>>>>>>>>> Convergence criterion not met.
 SCF Done:  E(RB+HF-LYP) =  -34352969.2846     A.U. after   65 cycles
             Convg  =    0.1924D-03             -V/T = *******
             S**2   =   0.0000
 Convergence failure -- run terminated.
 Error termination via Lnk1e in /sbin/loader.
 Job cpu time:  6 days 21 hours 27 minutes 49.9 seconds.
 File lengths (MBytes):  RWF=  539 Int=    0 D2E=    0 Chk=   30 Scr=    1

Any hints?  TIA

============================== Responses ==============================

From: Matt Challacombe T-12 <[email-masked]>

It's hard to know exactly what is happening, but you may
have instabilities from the "single point" cutoffs.  This can
cause small but serious electron leakages, as can a cheasy
grid.  I think you should tighten up your integrals, and try
again with a more verbose output.  Out of curriosity, what
kind of system is it?

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

From: Han Zuilhof <[email-masked]>

This frequently has happened to me as well, and a good work-around is
usually the following:
go from HF/3-21G* to B3LYP/3-21G* (just single point), then to
B3LYP/6-31G(d,p) (just single point again), and only finally to
B3LYP/6-31+G(d,p).

Use
geom=check & guess=read
in every step, and the same checkpoint-file.  I.e. a job like (for a
neutral singlet system):

%chk=marc
%mem=6000000
#p b3lyp/3-21G* geom=check guess=read

marc's structure

0,1

--link1--
%chk=marc
%mem=6000000
#p b3lyp/6-31G(d,p) geom=check guess=read

marc's structure

0,1

--link1--
%chk=marc
%mem=6000000
#p b3lyp/6-31+G(d,p) geom=check guess=read scf=tight

marc's structure

0,1

An approach like this worked in nearly all of my b3lyp/6-31+G(d,p) jobs!

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

From: Serguei Patchkovskii <[email-masked]>

Try to run it with SCF=(NOVARACC,NOINCFOCK,TIGHT), and do NOT
use GUESS=READ. Incidentally, your "very fast machine" is either
not very fast indeed, or heavily loaded. A 90MHz R8k PowerChallenge
(which is reasonably fast, but not the hottest machine out there
by a very long shot) can do single-point B3LYP/6-311G** calculation 
on a molecule with *120* carbon atoms in under a week.

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

From: Georg Schreckenbach <[email-masked]>

well, I guess everybody who tried to do metal systems with Gaussian has
seen that kind of convergence failure, at least I find it all the time. It
is, indeed, pretty annoying, and for such a large calculation in
particular. Here are a few things that may help.

1) The reason DFT is generally harder to converge than HF is that the
HOMO-LUMO gaps are larger (too large) in HF. Further, larger basis sets
with diffuse functions are often harder to converge than smaller ones. In
fact, this convergence failure is not all that surprising from my
experience.
2) Typically, I do, for metal systems that are hard to converge, a stepwise
procedure of "building the guess". It goes somehow like this.
a) Single point HF calculation at a small basis set level. You have done
that already.
b) Feed in the resulting wavefunction (guess=read in G94) into a single
point DFT (e.g., B3LYP) calculation at the same basis set level.
c) Use this result as input for a somewhat higher basis set (e.g., 6-31G or
6-31G*) ...
d) ... continue like that with successive single point calculations until
you reach the desired level of basis sets.
3) Further, I do, at a given level of basis set, two single point
calculations each: first with level shifting (e.g., Vshift=500), and second
without, again feeding the results of one calculation into the next one.

Whether all of this is necessary -- or feasible! -- for your system, I
don't know. It works, however, pretty well for a lot of different
molecules.

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

From: "Steven . Trohalaki MLPJ LAI" <[email-masked]>

marc - i had a similar problem but my molecule had a Ni atom.  still,
the solution to your problem may be the same.   doug fox said that
the problem was a poor initial guess and suggested that i use the
wavefunction from an HF run.  i'll append his original reply to me below.

steve t

   One problem relates to the fact that Gaussian needs a better initial guess
generator for metals.  This is in the works but will not be in G94.  The
other problem can be independent of initial guess, there are a large number
of potential low lying states for metal systems.  If you look at the SCF the
energy is oscillating wildly which indicates to me that simply letting it
iterate will not solve things and a better initial guess is required.

   First, move down to a single point with a smaller basis.  STO-3G or
more likely 3-21G.   LANL1MB is an ECP option.  These have fewer degrees
of freedom and so even a poor guess will converge.  Use POP=REG on
this so you can see the valence region and you can look for either
virtuals with negative energies or positive occupied, often signs that
GUESS=(READ,ALTER) might shift to a lower energy.  STABLE=OPT can be
helpful but this requires a fully converged SCF to begin with and even
with the smaller basis you might not get this.  Then use GUESS=READ
with the bigger basis set.

   Second, you might try using HF to generate the initial guess.  The
definition of the valence space is often quite different and cases where
SCF converges DFT does not.  Then use B3LYP GUESS=READ SCF=TIGHT and
see if you can get full convergence before launching the optimization.
You might also try SCF=QC at this point since the initial guess is fairly
close.

   Third, if you still have problems try running the dication.  This
overbinds the remaining electrons and can give you a well defined set
of core orbitals and with POP=REG followed by GUESS=(READ,ALTER) you
can pick the most likely orbital for the remaining two electrons.

  There is no magic formula and if there is no dominant single configuration
you might continue to have problems.  You might need to consider a
UHF singlet using GUESS=(READ,MIX) to break the spin symmetry which is
what STABLE=OPT may end up doing in a more formally correct form.

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

From: Krzysztof Radacki <[email-masked]>

First - it's error termineted because in 65 steps it couldn't find minima.
Second - I would presume  that difference between 631+** and 321* is big so
         you have not a  very nice starting point for your b3lyp.
Third - 631+** it's REALY high level. I did some opt. with 7 atoms from 
        second period (3* methyl group) on 631++** and it takes TIME.
Fourth - fast computer is relative ;) for me at home double-PII 300MHz 
         (ca. 300-400MFlops) would be very fast (I haven't it) on uni 
         Challenge M with R8000 processor is slow.  :(

Practicaly. I would sugest some medium level of complication. I never made 
321* with b3lyp but maybe it can help (or 631* if not)

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

From: Timothy Dransfield <[email-masked]>

I'm afraid I don't have a very good answer for *why* this works, and it's
not a particularly satisfying solution, but you can use SCF=QC for those
points where convergence is a problem.  I see this at long range, or in
extended bonds, myself.  Alas, SCF=QC makes the job take an order of
magnitude longer, so I'm hoping someone provides you with a better
workaround.  Please let me know if they do!!

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

From: "Oscar N. Ventura" <[email-masked]>

I have found the same problem in the study of transition metal complexes.
My workaround has been to use a smaller basis set. I employ
normally a B3LYP/3-21G* calculation after the HF one. I later perform the
B3LYP/6-311+G(2d,p) calculation or similar employing the GUESS=CHECK
option, so that the calculation starts from the converged B3LYP/3-21G*
density. The only caveat is that sometimes I have found that the
B3LYP/3-21G* calculation can converge to a different state of that I am
interested in. Just check the output.

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

From: Peter Sturm <[email-masked]>

Gaussian gives some hints about failed SCF convergence 
in its "Gaussian94 User's Reference" in Chapter 4, p. 162  under
"Problem Convergence Cases"

If you haven't this Manual, just look at

http://www.gaussian.com/g94.htm
 under "Using Gaussian 94 Efficiently" - "SCF Energies and Gradients"

or directly at

http://www.gaussian.com/00000029.htm

There you can find the right part of the online manual. The entry
"PROBLEM CONVERGENCE CASES" is identical to the text of the User's
Reference.

Perhaps you will find there a solution to your SCF convergence problem,
but trying out can be very time expensive!

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

From: Frank Jensen <[email-masked]>

        Try adding some level shift (vshift=xxx)

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

From: German Sastre Navarro <[email-masked]>

  it seems the initial guess from your previous calc is really far from
  being good. you can either try to:

  1./ increase the maximum scf cycles with: SCF+(MAXCYC=100)
      (i guess your default is 65)

  2./ try to re-run the calc without reading the scf from the previous calc

  I was anyway very surprised for this to take 6 days. May be you can speed
  it up with a bit more physycal memory.

  You can also try to run the calc with a less severe convergence criteria
  such as: 
        SCF=(conver=5)  sets conversion to 10^-5
        SCF=sleazy      reduce cutoffs for integrals
  otherwise you may end up with another 6-days run

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

From: Tapas Kar <[email-masked]>

It is sometime worth to use #P in the 1st line of com file to get more 
info on energies in each step. This gives an idea about the oscillating 
pattern of energies. Normally we use vshift if energy is oscillating. Try 
with more scf cycles, IOP(5/7=n), n is the no. of cycles you want. SCF=QC 
is another option but it takes huge CPU time.
You can restart the SCF from checkpoint file with vshift and increased 
number of scf steps. ( some cases i used 150 or even 
more number of cycles )

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

From: Thomas Dargel <[email-masked]>

in my experience, I would say Tapas is partially not right:

>You can restart the SCF from checkpoint file with vshift and increased
>number of scf steps. ( some cases i used 150 or even
>more number of cycles )

You could only restart and use the new produced and unconverged wafe
function if you plug iop(5/13=1) in the route card. This forces G94
to handle the unconverged SCF calculation not as an error, and writes
the gues back to checkpoint.
So, you can do several singlepoint calculations, which read in the
guess from before. This works really good for me if I have had a problem
with convergence in SCF for molecules containing transition metals.

E.g.:

%mem=...
%chk=/home/marc/dG-C8-FAAF_01-711_75.49.chk
#p B3LYP 6-31+G(d,p) SP scf=(nosp) scfcyc=25

Title

0 1
geometry

--link1--
%mem=...
%chk=/home/marc/dG-C8-FAAF_01-711_75.49.chk
#p B3LYP 6-31+G(d,p) SP scf=(nosp) scfcyc=25
geom=check guess=check

Title

0 1

--link1--
    .
    .
    .
    .

scfcyc=25, because I think the SCF also converges better using the
re-calculated integrals. It is better than running one job with
200 or more SCF cycles.

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

From: Bill Glauser <[email-masked]>

We came across your recent posting about the convergence problem encountered
while attempting a single point DFT calculation.

The Jaguar ab initio/DFT program that we develop and market is significantly
faster and more robust than other programs - especially for large systems.
It therefore occurred to us that Jaguar might be able to handle this case.

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

------------------------------------------------------------------------
 Marc C. Nicklaus                        National Institutes of Health
 E-mail: [email-masked]               Bldg 37, Rm 5B29
 Phone:  (301) 402-3111                  37 Convent Dr, MSC 4255       
 Fax:    (301) 496-5839                  BETHESDA, MD 20892-4255    USA 
      http://rex.nci.nih.gov/RESEARCH/basic/medchem/mcnbio.htm
    Laboratory of Medicinal Chemistry, National Cancer Institute,  &
  Center for Molecular Modeling, Ctr. for Information Technology, NIH
------------------------------------------------------------------------