Part of LOM

This blog is part of web activities of the Laboratory of Organic Materials (LOM) of the Institute of Solid State Physics of the University of Latvia.

Wednesday, December 28, 2016

Optimization with CalcAll not converged

Usually the scripts I use simply go to the next job step with Opt=CalcAll. Now, if THAT hasn't converged, it seems to be a serious trouble. I am far too lazy to change the initial geometry (although this is quite probable to help), so I try the following possibilities:
  • Opt=(CalcAll,NoRaman,MaxStep=2) (size of initial step taken by optimization algorithm; size of further steps is determined by the program as usual)
  • Opt=(CalcAll,NoRaman,RFO) (alternative algorithm)
  • Opt=(CalcAll,NoRaman,Cartesian) (optimization in Cartesian coordinates instead of redundant internals – works VERY rarely!)
  • Integral=UltraFineGrid (for meta functionals might help; also, makes sense because by default in G09 and G16 the grid used for CPHF/CPKS, which in involved in calculating the harmonic frequencies, is two steps behind the grid used for the density proper). If You do not want to repeat all other calculations with the better grid, this new one can be a better starting structure for calculation with G09 default FineGrid. (for CPKS, it's the Coarse). This problem is eliminated in G16, where the UltraFineGrid is the default (with SG1 for CPKS).
  • Also, Acc2E=11 if the CPKS problems are the culprit.
  • Sometimes switching from one algorithm to another in subsequent calculations may help.
  • With no other options, one have to use input file from different calculation method, different solvent or different ionization state.

First three options are easily implementable into Your scripts.

SCF Convergence failure. YQC and its friends

Lots of this already on the Internet, but for additional reference here is what I use in such cases, in the usual sequence:
  1. SCF=(YQC,IntRep,MaxConventionalCyc=40)
  2. SCF=(Fermi,IntRep,YQC,MaxConventionalCyc=40,NoVarAcc)
  3. SCF=(Fermi,IntRep,NoVarAcc)
  4. IOp(5/22=20) (pure EDIIS; by default, EDIIS is applied after CDIIS/EDIIS error becomes too big)
  5. Combinations of IOp(5/22=20) with all the previous variants.
IntRep is for symmetrical cases and usually shouldn't harm. If your molecule really has a symmetry other than C1, maybe Symmetry=Loose can help (usually input coordinates are a little distorted). NoVarAcc prevents too loose integral accuracy in early cycles of the SCF procedure (which usually speeds up the calculation, but may stab it in the back sometimes, too).

Note for code users:
Not always SCF failure is manifested by the

>>>>>>>>> Convergence criterion not met.

message. If the procedure failed at steepest descent stage, the following is printed:

No lower point found -- run aborted.

NB: in this case, there is no

Convergence failure -- run terminated.

message either.

Discussion of YQC output

YQC is an ingenious algorithm implemented in Gaussian. Citing the Release notes for Rev. D.01,
It does steepest descent and then scaled steepest descent as in SCF=QC, but then switches to regular SCF instead of quadratic convergence, using the quadratic algorithm only if the regular SCF fails to converge.
MaxConventionalCyc option set the number of cycles allowed for the direct SCF; usually if the calculation has not converged after 40 cycles, there is need for algorithm change. I do not remember where did I get the number, but most probably it is the ccl.net mailing list (amazing thing, check it out!).

One additional note for YQC: as the program exits one of the algorithms mentioned, it prints the SCF Done: line with corresponding energy and convergence. Usually the final energy appears the second SCF Done: line of YQC; this line can be distinguished by that the energy values are in capitals: "A.U." instead of "a.u.". Some people I know are used to capital "A.U." for finding the "right" energy for each optimization step, and usually it is completely OK. For example,

 SCF Done:  E(UCAM-B3LYP) =  -1440.14404401     a.u. after    7 cycles
 S**2 before annihilation     0.0000,   after     0.0000
 SCF Done:  E(UCAM-B3LYP) =  -1440.14408203     A.U. after   65 cycles
 S**2 before annihilation     0.0000,   after     0.0000
 QCSCF skips out because SCF is already converged.

But You should never use this in Your scripts. In rare cases, there might appear three subsequent SCF Done: lines, like here, for example:

 SCF Done:  E(UCAM-B3LYP) =  -1440.14404401     a.u. after    7 cycles
 S**2 before annihilation     0.0000,   after     0.0000
 SCF Done:  E(UCAM-B3LYP) =  -1440.14408203     A.U. after   65 cycles
 S**2 before annihilation     0.0000,   after     0.0000
 SCF Done:  E(UCAM-B3LYP) =  -1440.14408203     a.u. after    1 cycles
 S**2 before annihilation     0.0000,   after     0.0000

As You can check by running some simple calculations (e.g., CH4 with HF/6-31G(d)), in fact, lowercase "a.u." is the output of Link 508 which does both [scaled] steepest descent and quadratic (Newton–Raphson) SCF. At the same time, "A.U." is the output of Link 502, which does the direct SCF algorithm. Also, there are other output differences for the two Links:

Link 502:
SCF Done:  E(UCAM-B3LYP) =  -2125.66438397     A.U. after  129 cycles
            NFock=128  Conv=0.18D-06     -V/T= 2.0090
<Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.5000 <S**2>= 0.8613 S= 0.5542


Link 508:
SCF Done:  E(UCAM-B3LYP) =  -2125.66438397     a.u. after    1 cycles
            Convg  =    0.3116D-05                     1 Fock formations.
              S**2 =  0.8613                  -V/T =  2.0090
<Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.5000 <S**2>= 0.8613 S= 0.5542


So "NFock=" is not an indicator either (ah, how convenient would it be).

The point is that, in Your scripts, You should rather search for SCF Done: until something that always goes immediately after the last appearance of that phrase. Something like Population analysis or GradGradGradGrad, or The wavefunction if we are processing a Stable job (then You shouldn't search for Population analysis, it gets repeated before the result on stability) For example, if using sed:

sed '1,/The wavefunction\|Population analys\|GradGradGradGrad/p' filename | tac | grep 'SCF Don'

This works only for the first SCF energy. The general solution is, for instance,

sed -n '/SCF Don/ {:bulka s/\n/<br>/g; N; /\n SCF Don/{D};/The wavefunction\|Population analys\|GradGradGrad/!b bulka;s/<br>.*//g;p}' filename

Here, sed searches for lines with "SCF Don", sets label at them (label is named bulka), then replaces all newlines with <br> and reads in the next line. Hence, only two lines will be in the buffer (the pattern space, PS) at once. Then, if it founds that the second line in the  PS starts with " SCF Don", it deletes the first one. Next, sed returns to the label set (the bulka) unless the new line contains some of mentioned "aborting patterns" (The wafefunction, Population analys or GradGradGrad). After one of these has been reached, everything in the PS beginning with <br> is deleted, and the resulting lines (the desired final SCF energies) are printed.
Surely, this code is lame. But You should get the idea.

Fermi or level shifting (VShift)

Fermi option calls for Fermi temperature broadening during the initial iterations along with CDIIS and damping. As the broadening is removed in subsequent iterations, it should not affect the resulting energy, and it should still be comparable with results of calculations without the Fermi option. Damping is mixing of the old density matrix with the new one to reduce the amplitude of its oscillations during the early iterations of SCF procedure.[1]

VShift option shifts orbital energies (by default, it is 100 mHa, but You can set it with VShift=N where N is in millihartrees, mHa). This can aid, but the resulting SCF energy will change also; that's [probably] why, as the Gaussian 09 Reference says, the automatic archiving is disabled by this option. You should only use calculations with this option as an initial guess for subsequent 'normal' SCF job. I personally have never encountered a job that required this additional level-shifting headache, though.

[1] F. Neese, F. Wennmohs, et al. ORCA Manual, v. 3.0, p. 289 (even if You do not use alternative programs, it is highly educative to read their manuals).