Gaussian input section ordering: Gen + Generic + External iteration

Sometimes I need to calculate solvent-wise relaxed excited states with a basis set not built in Gaussian. Often when I am doing this I encounter an error: either the SCRF module is trying to interpret the basis as a solvent definition:

 Using the following non-standard input for PCM:
 C H N O 0
 ------------------------------------------------------------------------------
 The following wrong keyword(s) have been found:
 C
 
 
 
 0
 ------------------------------------------------------------------------------

  or the actual solvent definition is confusing the program: 

 Unrecognized atomic symbolNonEq

 So what is going on? If we look at the official Gaussian page for the Input section ordering [1]:

 

we see that there are two SCRF sections! The problem with our input file was that we assumed that the ExternalIteration keyword is everywhere while in reality it is not present in the first job. Typically this happens because of the copy-paste attitude, also when the input files are prepared by a [non-ideally written] script.

The way I know works is as follows:

1. For the first job with NonEq=Write we put the SCRF additional input after the Gen/GenECP additional input;
2. For the rest of jobs, we do the opposite.

I hope this helps :)

 

[1] https://gaussian.com/input/

ADMP error "IAtGet: Unknown atom type"

Now this is something that is not well documented by Gaussian, Inc. Description for the ADMP keyword states that most jobs will not require [additional input]. And yet if you have some additional input with lower priority than it would be for ADMP, the job dies and the following error appears:

IAtGet:  Unknown atom type

Considering that lower priority: look up the very beginning of the Gaussian NN Reference or go here (on the Gaussian homepage). For example,

• additional input for Opt=ModRedundant or CPHF=RdFreq is of higher priority than the one for ADMP, while the
• additional input for SCRF=Read or basis set specifications for the Gen keyword are of lower priority than the input for ADMP.

Actually, there is a hint on the mentioned pages of the Reference – the keyword field reads: ADMP and BOMD required input and ReadVelocity, ReadMWVelocity options, so we can actually figure out that there is some required input for ADMP, just like in the case of BOMD.

The solution is incredibly simple: you just have to add a blank line between the sections of additional input, the last one of higher priority and the first one with lower priority than the one of ADMP. For example, instead of

%Mem=120GB
%NProcShared=16
%Chk=/home/igors/darba/IN_dihedral_ADMP_DL1.chk

#T MaxDisk=250GB CAM-B3LYP/6-311G(d,p) ADMP=(MaxPoints=10000) Geom=(Crowd,ModRedundant) Guess=Mix scrf=(cpcm,solvent=Acetonitrile,Read)

::: IN17_469_CAM_6311_acetonitr_SASdrc_xplc_90deg IN17_469_CAM_6311_acetonitr_SASdrc_xplc_90deg GeometryOptimization Admp ADMP 1

0 1
   C     3.440990    2.083690   -0.808840
   C     3.774870    0.798710   -0.409860
   C     2.843130   -0.009000    0.245700

D 8 9 11 12 F

Surface=SAS Dis Rep Cav

you should have

%Mem=120GB
%NProcShared=16
%Chk=/home/igors/darba/IN_dihedral_ADMP_DL1.chk

#T MaxDisk=250GB CAM-B3LYP/6-311G(d,p) ADMP=(MaxPoints=10000) Geom=(Crowd,ModRedundant) Guess=Mix scrf=(cpcm,solvent=Acetonitrile,Read)

::: IN17_469_CAM_6311_acetonitr_SASdrc_xplc_90deg IN17_469_CAM_6311_acetonitr_SASdrc_xplc_90deg GeometryOptimization Admp ADMP 1

0 1
   C     3.440990    2.083690   -0.808840
   C     3.774870    0.798710   -0.409860
   C     2.843130   -0.009000    0.245700

D 8 9 11 12 F


Surface=SAS Dis Rep Cav

Yes, that's it. No more, no less. Everything goes on smoothly then.

Actually, this is in full analogy with BOMD, which requires the number of dissociation paths in the additional input section, but accepts blank line instead of zero.

P.S. Of course the structure in the input is not a real one :-)

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.

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.

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