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).

Thursday, November 17, 2016

That Dyson hand-dryer was not in the right place, from now on such images will appear in this page: http://cfilomquantum.blogspot.com/p/blog-page.html

Thursday, November 10, 2016

Differentiation step in Polar=Gamma

This is quite important note for those who wish to alter the defaults. 

Gaussian 09 Reference explains that the electric field step size is set by the option Step=N, by which the actual step size is set to 0.0001N atomic units (a. u.). This is, however, true only if you use double-numerical differentiation (done by Link 111). If You use numerical differentiation only once (like You do in Polar=Gamma calculations to obtain second hyperpolarizability γ), Link 106 is invoked instead, and it has different defaults. Namely, electric step field size is then 3 times smaller than it would be for Link 111, either by using the defaults or by reading user-specified values. So, when I had in the Route section
Polar=(DCSHG,Gamma,Step=3)
which, I thought, would do the calculation with the default field strengths of ±0.0003 a. u., I got instead
The following finite field(s) will be applied:  
An electric field of             -1.0000D-04  0.0000D+00  0.0000D+00
so -0.0001 a. u. Because Link 106 is invoked in Gamma calculations, and not the Link 111, for which finite-field parameters are specified by default!
In fact, this seems to me quite logical. In the "old times", finite-field (FF) and coupled-perturbed HF/KS (CPHF/CPKS) were used in standalone fashion as alternative methods for computing the hyperpolarizability β; then, of course, double-numerical differentiation was used, and maybe that is why the default "language" of parameter input is for the corresponding Link 111. Now, however, the second hyperpolarizability calculations are introduced into Gaussian, but only in a "hybrid" fashion of numerically differentiating β tensor. So, the Link 106 is now in use, but with different optimal parameters.

Sadly, these differences are listed only in the Gaussian IOp Reference. So, be aware of it and carry on Your research.

Friday, April 22, 2016

ntrex1 error of Gaussian


Error


ntrex1

as a separate line at the bottom of output file occured recently, as I was transferring an .GJF file from Windows™ onto Linux® and forgot to change path for all checkpoint files. So, the system could not find "D:\darba\something.chk", and the process crashed with the aforementioned error. Please check twice when You are changing Your checkpoints! :-)

Thursday, March 31, 2016

Gaussian performance on Windows and on Linux

As our Institute has licenses for Gaussian™ 09, Revision D.01 on both Windows™ and Linux®, it was interesting to compare performance on a single machine where both operating systems are installed. It has Intel® i3 2-core CPU (4 threads), 4 GB of RAM and Samsung ST500 hard drive – standard desktop PC, actually. Usually people test computational software on a single CPU core, but we ran both variants. Various calculation types were tested for some of our molecules.

At first glance, there is a big inconsistency in results: for single-core computations Gaussian 09 on Linux performs significantly better than on Windows™, but for multiple-core test the situation seems to be opposite. But if we divide "Job cpu time" values by the number of CPU treads for Linux values, the situation looks more logical... At least, these numbers are comparable.


Windows 7


Job type Calculation with single CPU core Calculation with 4 CPU cores Calculation with 4 CPU cores, norming to single core


Opt
Freq
Stable=Opt
Polar
Polar + SCRF
Total:

Hours Min. Sec.
2 30 40
6 56 41
1 23 33
1 28 37
1 42 0
14 1 31

Hours Min. Sec.
1 21 55
3 38 8
0 33 0
0 43 47
0 50 3
7 6 53

Hours Min. Sec.
0 20 28.75
0 54 32
0 8 15
0 10 56.75
0 12 30.75
1 46 43.25
Total: 50,491 sTotal: 25,613 sTotal: 6,403.25 s



Debian GNU/Linux 8.1


Job type Calculation with single CPU core Calculation with 4 CPU cores Calculation with 4 CPU cores, norming to single core


Opt
Freq
Stable=Opt
Polar
Polar + SCRF
Total:

Hours Min. Sec.
1 30 37.7
4 18 13.8
0 36 55.1
1 6 27.1
1 13 41
8 45 54,7

Hours Min. Sec.
3 33 27.3
11 29 34
1 18 14.8
2 21 4.9
2 34 58.3
21 17 19.3

Hours Min. Sec.
0 53 21.825
2 52 23.5
0 19 33.7
0 35 16.225
0 38 44.575
5 19 19.825
Total: 31,554.7 sTotal: 76,639.3 sTotal: 19,159.83 s


My personal conclusion is that on Unix™, Gaussian™ returns calculation time as if single CPU core was used on Unix™, but on Windows™, the actual computation time is returned. This was confirmed just by comparing time of creation and the last modify time for each file: on Linux this time span was far shorter than it would be if we sum up all "Job cpu times". Therefore, third section of the first table and second section in the second table are not corresponding to reality.

A practical conclusion is that on UnixGaussian™ runs faster than on Windows™. A technical conclusion is that we see clearly that not so much of computation can be parallelized, actually. Also, parallelization improves results much more on Windows than on Linux (but they are still worse).

We have also tried to use specific proprietary CPU firmware on Linux (Debian package i3fw). However, the results generally became slightly worse (for almost all jobs). We won't judge on it, because this "slightly" is really marginal difference. However, if You are free software freak, I think these findings will warm Your heart :)

Sunday, March 27, 2016

OrtVc1 failed #1.

Well, no much information about this on the internet. We have run into this when doing polarizability calculations on a supercell cut from a molecular semiconductor crystal.
Only information I have found is available here: http://www.somewhereville.com/?p=2175.

Disabling fast multipole method (by NoFMM keyword which is denoted as "obsolete" keyword in G09 Reference) helped us to resolve the issue. We, however, wrote to Gaussian Tech support and received further advice. The key point is that the problem arises if the system considered has some symmetry and it is used in the calculation; then FMM for electric field CPHF can fail. Fast multipole method is used only for large enough systems, and not only for CPHF; so, Gaussian' s Dr. Clemente advised me to disable it only for CPHF calculation by using IOp(10/63=1) keyword in the route section of Gaussian input file.

Gaussian technical support can be contacted through help@gaussian.com, usually they are very helpful.

Friday, March 25, 2016

RHF, UHF and Guess=Mix

It is straightforward to deduce that by default UHF (or UKS) singlet calculation produce identical alpha and beta orbitals, because there is no reason to break spin symmetry during the optimization (initially orbitals for both spins are the same, as follows from physical considerations). This was, however, new for me (thanks for explanation to Gaussian support team) when I had become suspicious about too good correspondence between RKS and UKS values, despite that I knew DFT is much less subject to spin contamination than HF.

Then, as it might happen if You are using some automation tool to generate input files, I have unintendedly generated some with RHF/... Guess=Mix. As in population analysis found in the output file there is only data on Alpha orbitals, I conclude that "R" or "U" in front of computation method is preferred over the Guess type by the program when it is determining the type of calculation.

UPDATE: First hyperpolarizability values are slightly changing for water and ozone calculations. RHF Guess=Mix seems to be closer to UHF than to RHF, although only alpha orbitals are shown in population analysis. Nevertheless, this difference is marginal. I will check later for some larger molecule.

Some default values in Gaussian 09

Probably these are considered "not interesting" by Gaussian, Inc., so not contained in G09 Reference (blue book). These and many others can be found in G09 IOps Reference (red book).

Integrals
2-Electron integral accuracy: 10-10
(i.e., Integral=(Acc2E=10) )

SCF Procedure
SCF convergence criterion: 10-8 , although for PBC 10-7
( SCF=(Conver=8) )
    Checked parameter is
  • RMS density for L502 (*DIIS)
  • RMS rotation gradient for L508 (Linear and Newton-Raphson quadratical convergence)
  • SQCDF for L506 (ROHF and GVB)
  • Energy for L510 (MCSCF).
RMS means "root-mean-square", not "Richard Matthew Stallman".

Polarizability
Numerical differentiation field step (which is quite important when You calculate Gamma) is by default 0.0003 a.u. (about 0.01542 V/Å). For double-numerical differentiation it is 0.001 a. u. (so 0.051 V/Å). For BOTH cases, this corresponds to keyword argument Step=10 (see this post).