TURBOMOLE input is keyworddirected. Keywords start with a ’$’, e.g. $title. Comments may be given after $dummy, or by a line starting with #; these lines are ignored by TURBOMOLE. Blank lines are also ignored. Keywords may be in any order unless stated otherwise below.
The sample inputs given below should help to give an idea how the keywords are to be used. They are sorted according to program. Complete control files are provided in Chapter 21. An alphabetical list of all keywords is given in the index.
The four keywords above are set by define, but are not necessary.
$statistics dscf
or
$statistics mpgrad
Only a statistics run will be performed to determine file space requirements as specified for dscf or mpgrad. On return the statistics option will be changed to $statistics off.
means current step. Keyword and data group (as e.g. dscf) is set by every program and removed on successful completion.
Keyword and data group (as e.g. relax) set by every program on successful completion.
General file crossreferences:
It is convenient not to include all input in the control file directly and to refer instead to other files providing the corresponding information. The above cross references are default settings from define; you may use other file names. define will create most of these files. Examples of these files are given below in the samples.
contains atom specification—type and location—and the bonds and internal
coordinates convenient for geometry optimizations.
specification of basis sets.
specification of effective core potentials.
MO vectors of SCF or DFT calculations for RHF or UHF runs.
keywords and data groups set by unrestricted dscf or ridft runs. Contain
natural MO vector and orbital occupation.
energies and gradients of all runs, e.g. for documentation in a geometry
optimizations.
approximate force constant for geometry optimizations.
The control file must end with this keyword:
General information defining the molecular system: nuclear coordinates, symmetry, basis functions, number of occupied MOs, etc. which are required by every module.
give title of run or project here.
Schönflies symbol of the point group. All point groups are supported with the
exception of NMR shielding and force constant calculations etc. which do not
work for groups with complex irreps (C_{3}, C_{3}h, T, etc). Use a lower symmetry
group in this case.
Example:
the basis set
and the auxiliary (fitting) basis for RIDFT calculations
the ECP if this is used.
The files basis, ecp and jbas must provide the necessary information under the labels specified in $atoms.
This data group specifies the number of Cartesian components of basis functions
(i.e. 5d and 7f in AOBasis, 6d and 10f in CAOBasis) for which the SCF calculation
should be performed. Possible values for char are AO (default) or CAO. If CAO is
used—which is not recommended—a core guess must be used instead of a Hückel
guess (see $scfmo).
Specification of MO occupation for RHF, e.g.
MO occupation of open shells and number of open shells. type=1 here means that
there is only a single open shell consisting e.g. of two MOs:
Roothaan parameters for the open shell, here a triplet case. define recognizes most
cases and suggests good Roothaan parameters.
For further information on ROHF calculations, see the sample input in Section 21.6 and the tables of Roothaan parameters in Section 6.3.
directs the program to carry out a UHF run,e.g.
The specification of MO occupation for UHF, $uhf overwrites closedshell occupation specification.
With the parameters in $redund_inp the generation of redundant internal coordinates can be modified. All entries have to be made in the control file before invoking the ired option. Important options are:
print parameter for debug output: The larger n is, the more output is printed
n ≥ 0,n ≤ 5 (default: 0)
method for generating and processing of redundant internal coordinates
n ≥  3,n ≤ 3,n ≠ 0 (default: 3)
Values for the metric option:
“Delocalized Coordinates”
The BmB^{t} matrix is diagonalized for the complete set of redundant
internal coordinates, matrix m is a unit matrix.
Delocalized Coordinates obtained with a modified matrix m, the values of m can be defined by user input (see below).
“Hybrid Coordinates”
Natural internal coordinates are defined as in the old iaut option. If
a cage remains, delocalized coordinates (as for n=1) are defined for
the cage.
Very simular to the n = 1 option, but for the remaining cage delocalized coordinates with modified matrix m are defined as for n = 3.
“Decoupled coordinates”
The redundant coordinates are divided into a sequence of blocks.
These are expected to have decreasing average force constants,
i.e. stretches, angle coordinates, torsions and “weak”
coordinates. The BB^{t} matrix is diagonalized for each block
separately after the columns of B were orthogonalized against the
columns of B of the the preceding blocks.
“Generalized natural coordinates”
Natural internal coordinates are defined first, for the remaining cage
decoupled coordinates are defined.
a positive real number, which is an approximate “force constant”, can be read in for
each type of coordinate (see below). The force constants are used for the definition
of the matrix m in BmB^{t}.
The matrix m is assumed to be a diagonal matrix. For each type of coordinate a different value for the force constants m_{ii} can be read in. Types of coordinates are:
bond stretch (default: 0.5)
inverse bond stretch (default: 0.5)
bond angle (default: 0.2)
Out of plane angle (default: 0.2)
dihedral or “torsional” angle (default: 0.2)
Special angle coordinate for collinear chains, bending of the chain a–b–c in the plane of b–c–d (default: 0.2)
bending of the chain a–b–c perpendicular to the plane of b–c–d
(default: 0.2)
stretch of a “weak” bond, i.e. the bond is assumed to have a very low force
constant, e.g. a “hydrogen bond” or a “van der Waals bond”
(default: 0.05)
inverse stretch of a weak bond (default: 0.05)
bond angle involving at least one weak bond (default: 0.02)
Out of plane angle for weak bonds (default: 0.02)
dihedral angle for weak bonds (default: 0.02)
linc coordinate for weak bonds (default: 0.02)
linp coordinate for weak bonds (default: 0.02)
One has to specify only the Cartesian coordinates (data group $coord) to start a uff run. The program uff requires the data groups $uff, $ufftopology, $uffgradient and $uffhessian. If these keywords do not exist in the control file the program will generate these data groups.
The data group $uff contains the parameters described below. The default values—in the control file—are:
The explanation of the variables are as follows:
number of max. optimization cycles (maxcycle=1: singlepoint calculation).
can have the values +1 or 1. If modus = 1 only the topology will be
calculated.
each nqeq cycle the partial charges will be calculated. If nqeq = 0, then the partial charges are calculated only in the first cycle, if the file ufftopology does not exist.
switch for the different types of force field terms:
bond terms will be calculated.
angle terms will be calculated.
torsion terms will be calculated.
inversion terms will be calculated.
non bonded van der Waals terms will be calculated.
non bonded electrostatic terms will be calculated.
convergence criteria for energy and gradient.
total charge of the molecule.
distance parameter to calculate the topology. If the distance between the atoms I and J is less than the sum of the covalent radii of the the atoms multiplied with dfac, then there is a bond between I and J.
if the norm of the gradient is greater than epssteep, a deepestdescentstep will be
done.
if the norm of the gradient is smaller than epssearch, no linesearch step will be
done after the Newton step.
max. displacement in a.u. for a coordinate in a relax step.
parameters of linesearch:
start value
increment
number of energy calculations
modification parameter for the eigenvalues of the Hessian (see below):
f(x) = x * (alpha + beta* exp(gamma* x)).
a switch for the transformation in the principal axis system.
a switch for the numerical Hessian.
a switch for a MD calculation.
cartesian coordinates of the atoms (default: $coord file=coord)
contains a list of the next neigbours of each atom (see Section 20.2.4).
Sometimes it is useful to enter the connectivity (in the input block nxtnei12 in
the file ufftopology) by hand (not always necessary; default: $ufftopology
file=ufftopology).
Beyond this uff reads the force field parameters for the atoms from the file parms.in. If this file exists in the directory from which one starts an uff calculation the program will use this file, if not the program reads the data from the file $TURBODIR/uff/parms.in. If one wants own atom types, one has to add these atoms types in the file parms.in. For each new atom type one has to specify the natural bond distance, the natural bond angle, the natural nonbond distance, the well depth of the LennardJones potential, the scaling factor ζ, the effective charge, torsional barriers invoking a pair of sp^{3} atoms, torsional barriers involving a pair of sp^{2} atoms, generalized Mulliken–Pauling electronegativities, the idem potentials, characteristic atomic size, lower bound of the partial charge, upper bound of the partial charge. Distances, energies and charges are in atomic units and angles are in rad.
contains the (updated) cartesian coordinates of the atoms (default: $coord
file=coord).
contains the full information of the topology of the molecule and the whole
force field terms (see below; default: $ufftopology file=ufftopology).
contains the accumulated cartesian analytical gradients (default:
$uffgradient file=uffgradient).
contains the cartesian analytical Hessian;
(default: $uffhessian file=uffhessian00).
The topology file ufftopology contains the blocks nxtnei12, nxtenei13, nxtnei14, connectivity, angle, torsion, inversion, nonbond and qpartial. It starts with $ufftopology and ends with $end. The first three blocks (nxtnei12, nxtnei13, nxtnei14) have the same form: they start with the atom number and the number of its neighbours, in the next line are the numbers of the neighbour atoms. Then the connectivityblock follows starting with the number of bond terms. Each line contains one bond term:

Here are I and J the number of the atoms, d the distance in a.u. and BO is the bond order.
The angle terms follow, starting with the number of the angle terms. In each line is one angle term:

Here are J,I and K the atoms number, where atom I is in the apex. “wtyp” is the angle type and has the following values:
linear case
trigonal planar case
quadratic planar case
octahedral case
all other cases.
θ is the angle value in degree. nr_{JI} and nr_{IK} are the number of the bonds between J and I and the bond between I and K. The hybridization of atom I determines “wtyp”.
Then the torsion terms follow, starting with the number of the torsion terms. Each line contains one torsion term:

Here are I,J,K and L the atom numbers. nr_{JK} is the number of the bond between J and K. “ttyp” is the torsion type:
J (sp^{3})–K (sp^{3})
like ttyp=1, but one or both atoms are in Group 16
J (sp^{2})–K (sp^{3}) or vice versa
like ttyp=2, but one or both atoms are in Group 16
like ttyp=2, but J or K is next a sp^{2} atom
J (sp^{2})–K (sp^{2})
all other cases.
ϕ is the value of the torsion angle in degree. θ_{IJK} is the angle value of (I  J  K) and θ_{JKL} is the cwone for J  K  L. The hybridizations of J and K determine “ttyp”.
The inversion terms follow starting with the number of inversion terms (e.g. the pyramidalisation of NH_{3}). In each line is one inversion term:

I,J,K and L are the atom numbers. Atom I is the central one. ityp1, ityp2, ityp3 are the types of the inversions:
atom I is C and atom L is O
like ityp=10, but L is any atom
I is P
I is As
I is Sb
I is Bi
all other cases.
ω_{1},ω_{2} and ω_{3} are the values of the inversion angles in degree.
The nonbond terms follow starting with the number of the nonbonded terms. In each line is one nonbond term:

Here I and J are the atom numbers, d the distance in a.u. Then the partial charges follow.
If the determination of the molecule connectivity failed, you can specify the block nxtnei12 in the file ufftopology. Then the program calculates the other blocks.
Based on the numbers of the next neighbours (block nxtnei12 in the file ufftopology) the program tries to determine the UFF type of an atom. The following rules are implemented: If the atom has three next neighbours and it is in the nitrogen group, then it has a hybridization three. If it is not in the nitrogen group, it has hybridization two. If the atom has four next neighbours and it is in the carbon group, it has hybridization three. If it is not in the carbon group, it becomes hybridization four. If the number of next neighbours is six, then it gets the hybridization six.
Since the smallest eigenvalues λ_{i} of the Hessian has the greatest influence on the convergence of the geometry optimization, one can shift these values with

and calculates a new Hessian with these modified eigenvalues.
Module WOELFLING reads options from data group
The below values of the options are default values with the following meaning:
Number of interpolated structures for optimization.
Number of input structures provided by user.
Align input structures by translation/rotation 0=yes, 1=no.
Maximum number of iterations.
Threshold for accuracy of LSTinterpolation.
Threshold for mean of norms of projected gradients.
Use standard optimization from initial LSTpath (method q) or grow reaction path (method qg).
Furthermore,
counts the number of completed iteration (no option).
Convergency criterion for the root mean square of the density matrix. If you
want to calculate an analytical MP2 gradient (program mpgrad) real should be
1.d7 or less.
Listing of all possible subkeywords for $dft (crossreferences are given).
The user normally has to choose only the functional and the grid size, see below. All other parameters have proven defaults.
Specification of the functional, default is BP86, printed as functional
bp. For all possible—and useful—functionals, please refer to page 653
and for definition of the functionals the section 6.2 on page 238.
Example (default input):
Specification of the spherical grid (see section 20.2.6 on page 653). Default is
gridsize m3.
Example:
—not recommended for use—
Specification of the mapping of the radial grid.
Possible values for gridtype are 1,…,6, but gridtype 4 to 6 is only for the use with functional lhf (see page 658). For the definition of gridtype 1 – 3, please refer to Eq. (16), (17) and, (19) in Ref. [192].
Example (default value):
Flag for debugging. debug 0 means no debug output (default), debug 1 means
some output, debug 2 means a lot more output. Be careful!
Specification of the sharpness of the partition function as proposed by
Becke [193], default is nkk 3. The usage of nkk makes sense only in the range
1 ≤ nkk ≤ 6.
Example:
—not recommended for use—
Only for userspecified Lobatto grids, i.e. gridsize 9, ntheta specifies the
number of θ points and nphi specifies the number of ϕ points. For the fixed
Lobatto grid, i.e. gridtype 8, the default value is ntheta 25 and nphi
48.
When gridsize 9 is given, you have to specify both, ntheta and nphi (see below), otherwise the program will crash. The restriction for userdefined Lobatto grids is the number of grid points, which must not exceed 2000 grid points.
Example:
Original grids had not been carefully optimized for all atoms individually. This
has now been done, which let to changes of ξ for Rb and Cs only resulting in
minor improvements. If you have ongoing projects, which have been started
with the old grids, you should continue using them with the keyword
old_RbCs_xi.
Example:
Specifies the number of radial grid points. Default values depend on type of
atom and grid (see keyword gridsize). The formula for the radial gridsize is
given as,

ioffrad is atomdependent, the more shells of electrons, the larger ioffrad:
elements  ioffrad  elements  ioffrad  
for H–He  20  for K–Kr  40  
for Li–Ne  25  for Rb–Xe  45  
for Na–Ar  30  for Cs–Lw  50 
The radial grid size increases further for finer grids:
gridsize  1  2  3  4  5  6  7  8  9 
radsize  1  2  3  6  8  10  14  9  3 
If you want to converge results with respect to radial grid size, increase radsize in steps of 5, which is convenient (see equation above).
Serves to increase quadrature grids; this is recommended in case of very diffuse
wavefunctions. With the keyword diffuse grids are modified by
changing the linear scaling factor ξ of radial grid points and the radial
gridsize:
radsize radsize + incr
ξ ξ * scal
diffuse integer  1  2  3  4  5  6 
incr  1  2  3  4  5  6 
scal  1.5  2.0  2.8  4.0  5.0  6.0 
For information about the linear scaling parameter ξ, see Eq. (16)–(19) and Table 1 in Ref. [192].
In addition, the reduction of spherical grid points near nuclei is suppressed, i.e. fullshell on is set (see page 627).
Note: the keyword radsize integer overrules the setting of incr; for more information, see p. 622.
Recommendation: For diffuse cases use gridsize m4 (or larger) in combination with diffuse 2 and check the number of electrons; for more difficult cases use diffuse 4. In case of doubt, verify the calculation with a larger grid, i.e. gridsize 7.
The test suite example $TURBODIR/TURBOTEST/dscf/short/H3PO4.DSCF.DIFFUSE provides an example of usage; this also gives reasonable values for damping and orbitalshift to reach convergence in this and similar cases, see $scfdamp and $scforbitalshift (p. 638 and p. 644).
Example (Recommendation):
—for developers only—
Radial grid points have a linear scaling parameter ξ, see Eq.(16)–(19) and
Table 1 in Ref. [192]. With the following input,
one performs a numerical integration for the density and the exchange correlation term for ξ = 0.5,(0.01),2.0 for given MOs and functional. NOTE: only molecules with a single atom type can be used. The results serve to establish stable, optimal ξ values, see Figure 1 in Ref. [192]. Program stops after this testing.
Usage of the reference grid, which is a very fine grid with very tight
thresholds. The default values for the different variables are:
Please refer to the corresponding subkeywords for explanation.
If you want to use the reference grid, you have to skip the keyword gridsize, and type instead reference. Example:
Check if the selected grid is accurate enough for the employed basisset by
performing a numerical integration of the norm of all (occupied and virtual)
orbitals. Useful for LHF. 659.
Grid points are sorted into batches, which are then processed. This increases
efficency. This should be changed only by developers. Default is batchsize
100.
Standard grids have reduced number of spherical grid points near nuclei.
With the keyword fullshell this reduction is suppressed. Reference grid
(see keyword reference) always has full spherical grids with 1202
points. Should be used to checked the influence of spherical grid
reduction.
Example for the usage of fullshell:
—for developers only—
Values of real effects efficiency of the quadrature, default is symblock1 0.001
and symblock2 0.001, one can try higher or smaller values.
—not recommended for use—
Where xparameter (default) can be: sgrenze (8), fgrenze (10), qgrenze (12),
dgrenze (12) and, fcut (14). These parameters control neglect of near zeros of
various quantities. With xparameter integer one changes the default. integer
larger than defaults will increase the numerical accuracy. Tighter threshold are
set automatically with keyword $scfconv (see section 20.2.6 on page
638).
Includes the derivatives of quadrature weights to get more accurate results.
Default is that the derivatives of quadrature weights will be not considered,
see section 20.2.9 on page 685.
Grid points are ordered into batches of neighbouring points. This increases
efficiency, since now zeros can be skipped for entire batches. gridordering is
default for serial version, not for the parallel one. You cannot use weight
derivatives and gridordering together.
Example for switching off gridordering:
Specification of external electrostatic field(s). The specification may take place either
by Ex, Ey, Ez or by x, y, z, E. See also $fldopt.
Example:
Requests calculation of occupation numbers at a finite temperature T. For an
orbital with the energy ε_{i} the occupation number n_{i} ∈ is calculated
as
where μ is the Fermi level. The factor f = 4k∕ is chosen to yield the same slope at the Fermi level as the Fermi distribution.
Calculation of the fractional occupation numbers starts when the current HOMOLUMO gap drops below the value given by hlcrit (default: 0.1). The initial temperature given by tmstrt (default: 300K) is reduced at each SCF cycle by the factor tmfac (default: 1.0) until it reaches the value tmend (default: 300K). Note that the default values lead to occupation numbers calculated at a constant T = 300K. Current occupation numbers are frozen if the energy change drops below the value given by stop (default: 1⋅10^{3}). This prevents oscillations at the end of the SCF procedure.
Calculation of fractional occupation numbers often requires much higher damping and orbital shifting. Please adjust the values for $scfdamp and $scforbitalshift if you encounter convergence problems.
In UHF runs this option can be used to automatically locate the lowest spin state. In order to obtain integer occupation numbers tmend should be set to relatively low value, e.g. 50K.
Calculation of fractional occupation numbers should be used only for single point calculations. When used during structure optimizations it may lead to energy oscillations.
The optional nue value (number of unpaired electrons) can be used to force a certain multiplicity in case of an unrestricted calculation. nue=0 is for singlet, nue=1 for dublet, etc.
Perform firstorder SCFcalculation, i.e. perform only one SCFiteration with the
start MOs (which should be the orthogonalized MOs of two independent subsystems
as is explained in detail in Chapter 16).
Specification of options related with external electrostatic fields. The following
options are available:
Calculate numerical 1st derivative of SCF energy with respect to
electrostatic field (default: off), increment for numerical differentiation
is edelt (see below).
Calculate numerical 2nd derivative of SCF energy with respect to
electrostatic field (default: off), increment for numerical differentiation
is edelt.
Increment for numerical differentiation (default: 0.005).
Calculate SCF energy for nonzero external electrostatic
fields defined in $electrostatic field.
Calculate SCF energy for one external field definition and dump
disturbed MOs onto $scfmo. This enables to evaluate properties or
perform geometry optimizations in the presence of an external field.
Caution: don’t use the RI approximation for all these calculations since this will lead to nonnegligible errors!!
By using this option the twoelectron integrals are kept in RAM; integer specifies
how many megabytes should be allocated. If the integrals exceed the RAM allocated
the program reverts to the standard mode. Supports all methods which process
twoelectron integrals, i.e. SCF and DFT (including hybrid functionals); RHF and
UHF.
The following condition must be met:
and rhfshells 1 or 2. It is advisable to set $thize as small as possible (e.g. $thize
0.1d08) and to remove the keyword $scfdump.
Note: this keyword does not work for parallel runs.
If this keyword is set the energies and symmetry labels of all occupied MOs will be
dumped to this data group. This may be helpful to draw modiagrams.
If only has been set only the start MOs are dumped and the program
quits.
nirreps will hold the total number of displayed orbitals after the successful run.
If this keyword is present all occupied orbitals are dumped to standard output. Be
careful about this option as it can create huge output files in case of many basis
functions.
If this line is present, the dscf program is forced to output the MOs using the new
FORTRAN format format regardless of the formatoption in data group $scfmo.
Otherwise the input format will be used.
Example: $mo output format(3(2x,d15.8))
This data group will be written after an UHF calculation (together with $natural
orbital occupation) and contains the natural space orbitals (same syntax as $scfmo).
This data group will be written after an UHF calculation (together with $natural
orbitals) and contains the occupation of natural orbitals (syntax as any data
group related with orbital occupation information, e.g. $closed shells),
e.g.
Specification of position and magnitude of point charges to be included in the
Hamiltonian. Each point charge is defined in the format
with <x>, <y>, <z> being the coordinates and <q> its charge,e.g.
In addition the following optional arguments may be given:
distance threshold for discarding redundant point charges, default value
10^{6}.
if given, the selfenergy of the point charge array will will be included in
the energy and the gradient
switches off the check for redundant point charges and the default
symmetrization. This option can significantly speed up the point charge
treatment if many of them are involved  use only if the point charges are
well distributed and symmetry is C_{1}, e.g. when they stem from proper
MM runs
print all point charges in the output (default is to print the point charges
only if less than 100 charges given)
concerns the first SCF iteration cycle if start MOs from an EHT guess are
used.
The SCF iteration procedure requires control mechanisms to ensure (fast) convergence, in TURBOMOLE these are based on orbital energies ϵ_{i} of the preceeding iteration used for level shifting and damping (besides DIIS, see below). This feature cannot be used in the first iteration if EHT MOs are employed as start, since ϵ_{i} are not available. The keyword $prediag provides ’ϵ_{i} of the zeroth iteration’ by diagonalization of occ–occ and virt–virt part of the first Fock matrix, to allow for level shifting etc.. See $scfdiis below.
Try a dscf restart. The program will read the data group $restartd (which must
exist, also $scfmo has to exist!) and continue the calculation at the point where it
ended before. If the additional option twoint is appended, the program will read the
twoelectron integrals from the files specified in $scfintunit, so there will be almost
no loss of cputime.
All this information is normally provided by the previous dscf run if the keyword $scfdump (see there) was given.
Data provided by a previous dscf run that has been interrupted. This keyword is
created when $scfdump was given.
is set by define so usually no changes are necessary. The dimensions must
be greater or equal to those actually required, i.e. you can delete basis
functions and keep rundimensions. This keyword is not necessary for small
cases.
Example:
SCF convergency criterion will be 10^{integer} for the energy. Gradients will only be
evaluated if integer > 6.
Damping parameters for SCF iterations in order to reduce oscillations. The old
Fockoperator is added to the current one with weight 0.5 as start; if convergence is
good, this weight is then reduced by the step 0.05 in each successive iteration
until the minimum of 0.1 is reached. (These are the default settings of
define for closedshell RHF). DSCF automatically tries to adjust the weight
to optimize convergence but in difficult cases it is recommended to start
with a large weight, e.g. 1.5, and to set the minimum to a larger value,
e.g. 0.5.
Flags for debugging purposes. Following options are available:
Output level concerning molecular orbitals. integer=0 (default) means
minimal output, >1 will output all start MOs and all MOs in each
iteration.
Output level concerning difference density matrices.
integer > 0 will dump a lot of information—be careful!
Direct SCF procedures build the Fock matrix by exploiting information from
previous iterations for better efficiency. With this keyword information
from the last integer iterations will be used. This feature is switched on
with the default value 20, even if the keyword is absent. The user may
reduce the number of iterations employed to smaller values (e.g. 10) in
cases were numerical stability could become an issue. With the value 0 this
feature is switched off; the Fock matrix is constructed from scratch in each
iteration.
Control block for convergence acceleration via Pulay’s DIIS
^{*}.
Options are:
specifies the kind of error vector to be used (two different kind of DIIS algorithms)
uses (FDS  SDF) as error vector.
no DIIS
use S^{1∕2}FDS^{1∕2}  transposed
Further suboptions:
maximum number of iterations used for extrapolation.
debug level (default: 0)
print applied DIIS coefficients
print DIIS matrix and eigenvalues, too
scaling factor in DIIS procedure: qscal > 1 implies some damping,
qscal = 1.0: straight DIIS.
directs the reduction of qscal to qscal = 1.0 (no damping in DIIS), done
if errvec < thrd.
Defaults for $prediag (see above) and $scfdiis
errvec=FDSSDF, maxiter=5, qscal=1.2, thrd=0.0, this implies DIIS damping in all iterations, prediag is switched of.
Recommended:
errvec=sFDs leads to the following defaults:
qscal=1.2, for SCF runs: maxiter=6 and thrd=0.3, prediag is off; for DFT runs:
maxiter=5 and thrd=0.1 prediag is on. If you want to switch off prediag put
$prediag none.
Dump SCF restart information onto data group $restartd and dump SCF MOs in
each iteration onto $scfmo (scfdump = iter). Additionally, a data block
$scfiterinfo will be dumped containing accumulated SCF total, one and
twoelectron energies of all previous SCF iterations. Information that will allow you
to perform a restart if your calculation aborts will be dumped on data group
$restartd (see also $restart).
Disc space specification for twoelectron integrals. The following suboptions are
available (and necessary):
Fortran unit number for this file. Unit numbers 30,31,… are
recommended.
Filespace in megabytes for this file. size=0 leads to a fully direct run.
size is set by a statistics run, see $statistics. DSCF switches to direct
mode if the file space is exhausted.
Filename. This may also be a complete path name, if you want to store
the integrals in a special directory. Make sure the file is local, otherwise
integrals are transmitted over the network.
Thus your data group $scfintunit may look like this:
Maximum number of SCF iterations (default: 30).
Input/output data group for SCF MOs. You can specify
To perform a calculation without a start vector (i.e. use a core
Hamiltonian guess).
The file where the MOs are written on output (default: mos).
These two options can also be used for $uhfmo_alpha and $uhfmo_beta to use a core guess and write the molecular orbitals to file.
After running define or a TURBOMOLE calculation additional options may appear specifying the origin of the MOs:
These MOs were obtained by projection form another basis set. They
should not be used for wavefunction analysis.
The MOs are converged SCF MOs , the convergence criterion applied
was 10^{integer}
The MOs are unconverged SCF MOs which were written on this data
group after iteration integer. The latter three options are mutually
exclusive.
This specifies the FORTRAN format specification which was used for MO
output. The standard format is (4d20.14). (See data group $mo output
format.)
Example:
Your data group $scfmo could look like this after a successful TURBOMOLE run
:
To assist convergence, either the energies of unoccupied MOs can be shifted to
higher energies or, in openshell cases, the energies of closedshell MOs to lower
energies. In general a large shift may help to get better convergence.
Options are:
Automatic virtual shell shift switched off.
Automatic virtual shell shift switched on; the energies of virtual orbitals
will be shifted if the HOMOLUMO gap drops below real such that a gap
of real is sustained. This is the default setting if the keyword is missing
with real=0.1.
Option for openshell cases. Closed shells are shifted to lower energies by
real. The default shift value is closedshell=0.4.
Note: Normally this will disable the automatic shift of energies of virtual
orbitals! To override this, you should append an exclamation mark to
the ’automatic’ switch, i.e. specify ’automatic! real’.
Set shifts for special occupied MOs. To use this option, start the line with the
symmetry label and the list of MOs within this symmetry and append the
desired shift in brackets as in the following example:
Integral evaluation threshold. Integrals smaller than real will not be evaluated. Note
that this threshold may affect accuracy and the convergence properties if it is chosen
too large. If $scftol is absent, a default value will be taken obtained from $scfconv
by real = (#bf = number of basisfunctions).
The scratch files allocated by dscf can be placed anywhere in your file systems
instead of the working directory by referencing their pathnames in this data group.
All possible scratch files are listed in the following example:
The following options are allowed
Do not perform integrals statistics
see KORA
see POLLY
see PARALLEL PROCESSING
Options kora, dscf parallel, grad, mpgrad, polly will be described in the
related chapters.
If $statistics dscf has been given integral prescreening will be performed (which
is an n^{4}step and may therefore be timeconsuming) and a table of the number of
stored integrals as a function of the two parameters $thize and $thime will be
dumped. Afterwards, the filespace needed for the current combination of $thize and
$thime will be written to the data group ($scfintunit) and $statistics dscf will
be replaced by $statistics off.
Integral storage parameter, which is related to the time needed to calculate the
integral. The larger integer the less integrals will be stored. The default value is
integer = 5. (see also $thize, $statistics)
Integral storage parameter, that determines, together with $thime, the number of
integrals stored on disc. Only integrals larger than real will be stored. The default
value is real = 0.100E04.
Specification of MO occupation for RHF, e.g.
MO occupation of open shells and number of open shells. ’type=1’ here means that
there is only a single open shell consisting e.g. of two MOs:
This data group is necessary for ROHF calculations with more than one open shell.
Example:
For ROHFcalculations with only one open shell the Roothaan
parameters^{†}
a and b have to be specified within this data group (see also $rohf). Example:
For further information on ROHF calculations (e.g. with more than one open shell), see the sample input in Section 21.6 and the tables of Roothaan parameters in Section 6.3.
Note that this keyword toggles the ROHF mode also for more than one open shell. If it is not given, the openshell electrons are simply ignored.
these two data groups specify the occupation of alpha and beta spin UHF MOs
(syntax as any data group related with orbital occupation information, e.g. $closed
shells)
Example:
directs the program to carry out a UHF run. $uhf overwrites closedshell occupation
specification.
These two data groups contain the UHF MO vectors for alpha and beta spin
respectively (same syntax as $scfmo).
for DFT calculations one has to specify the functional and the grid (for the quadrature of the exchange correlation part). The settings above are default: both lines can be left out if the BP86 functional and grid m3 are required. Other useful functionals supported are:
(equivalent to the Gaussian98 keyword B3LYP with VWNIII)
(equivalent to the Gaussian98 keyword SVWN with VWNIII)
Possible grids are 1–5 and m3–m5 where grid 1 is coarse (least accurate) and 5 most dense. We recommend however the use of socalled multiple grids m3–m5: SCF iterations with grid 1–3, final energy and gradient with grid 3–5. Usually m3 is fine: for large or delicate systems, try m4. For a reference calculation with a very fine grid and very tight thresholds use ’reference’ as grid specification instead of ’gridsize xy’.
Note: the functionals b3lyp_Gaussian and svwn_Gaussian are made available only for comparability with Gaussian. The functional VWNIII is much less well founded than VWN5 and the TURBOMOLE team does not recommend the use of VWNIII.
Dscf does not run with the keyword $rij: you must call the RI modules Ridft and Rdgrad for energy and gradient calculations. However, it does run with the keyword $rik, but it will ignore all RI settings and do a conventional nonRI Hartree–Fock or DFT calculation.
Enforces an RIJ calculation if module ridft is used, can be used for
HartreeFock as well as for DFT calculations with pure or hybrid functionals.
Obsolete keyword  use $rij instead!
Enforces a RIJK calculation if module ridft is used, can be used for
HartreeFock as well as for DFT calculations with pure or hybrid functionals.
Choose the memory core available (in megabyte) for special arrays in the RI
calculation (the less memory you give the more integrals are treated directly,
i.e. recomputed on the fly in every iteration)
Cross reference for the file specifying the auxiliary basis as referenced in
$atoms. We strongly recommend using auxbasis sets optimized for the
respective MO basis sets, e.g. use SVP (or TZVP) for the basis and the
corresponding auxbasis as provided by define (default: file=auxbasis).
Calculation of atomic charges according to the s partial wave and atomic
dipole moments according to the p partial wave as resulting from the auxbasis
representation of the density
If the keyword $rik is found in the control file, ridft performs a Hartree–Fock–SCF calculation using the RIapproximation for both Coulomb and HFexchange (efficient for large basis sets). For this purpose needed (apart from $ricore):
Cross reference for the file specifying the JKauxiliary basis as referenced in
$atoms. This group is created by the rijk menu in define.
MultipoleAcceleratedResolutionofIdentityJ. This method partitions the Coulomb interactions in the near and farfield parts. The calculation of the farfield part is performed by application of the multipole expansions and the nearfield part is evaluated employing the RIJ approximation. It speeds up calculation of the Coulomb term for large systems. It can only be used with the ridft module and requires setting of the $ridft keyword.
The following options are available:
specifies precision parameter for the multipole expansions. Lowprecision MARIJ calculations require 1 ⋅ 10^{6}, which is the default. For higher precision calculations it should be set to 1 ⋅ 10^{8}–1 ⋅ 10^{9}.
maximum lmoment of multipole expansions. It should be set to a value equal at least twice the maximum angular momentum of basis functions. Default value is 10 and it should probably never be set higher than 18.
Threshold for moment summation. For highly accurate calculations it should be set to 1 ⋅ 10^{24}.
number of bins per atom for partitioning of electron densities. Default value is 8 and hardly ever needs to be changed.
minimum separation between bins. Only bins separated more than the sum of their extents plus wsindex are considered as farfield. Default is 0.0 and should be changed only by the experts.
maximum extent for charge distributions of partitioned densities. Extents with values larger then this are set to extmax. Hardly ever needs to be changed.
If the keyword $senex is found in the control file, ridft performs a Hartree–Fock– SCF calculation using the seminumerical approximation for HFexchange. Standard dftgrids can be used for the numerical integration. Smaller grids (1,0) and the corresponding mgrids (m1,m2) have been defined as well. Grids of at least size m3 are recommended for heavy atoms. The gridsize can be modified just like in dftcalculations. The keyword $dsenex activates seminumerical gradient calculations. An example using the default grid for SCF (m1) and grid 5 for gradients (default grid: 3) looks like this:
Use the Localized Hartree–Fock (LHF) method to obtain an effective ExactExchange Kohn–Sham potential (module dscf). The LHF method is a serial implementation for spinrestricted closedshell and spinunrestricted ground states.
With the LHF potential Rydberg series of virtual orbitals can be obtained. To that end, diffuse orbital basis sets have to be used and special grids are required.
gridtype 4 is the most diffuse with special radial scaling; gridtype 5 is for very good Rydberg orbitals; gridtype 6 (default in Lhfprep) is the least diffuse, only for the first Rydberg orbitals.
Only gridsize 3–5 can be used, no multiple grids.
Use testinteg to check if the selected grid is accurate enough for the employed basisset.
Otherwise the LHF functional can be selected in define: in this case default options are used.
Options for the LHF potential can be specified as follows ( see also lhfprep help)
calculation of the KLI exchange potential. By default the LHF exchange
potential is computed (offdiag on).
the Slater potential is calculated numerically everywhere: this is more
accurate but much more expensive. When ECP are used, turn on this
option.
leads to accurate results only for firstrow elements or if an uncontracted
basis set or a basis set with special additional contractions is used: in
other cases numericalslater on has to be used (this is default).
for asymptotic treatment there are three options:
No asymptotictreatment and no use of the numerical Slater. The
total exchange potential is just replaced by 1∕r in the asymptotic
region. This method is the fastest one but can be used only for the
densitymatrix convergence or if Rydberg virtual orbitals are of no
interest.
Full asymptotictreatment and use of the numerical Slater in the
near asymptoticregion.
Automatic switching on (off) to the special asymptotic treatment
if the differential densitymatrix rms is below (above) 1.d3. This
is the default.
the converged Slater and correction potentials for all grid points are saved in
the files slater.pot and corrct.pot, respectively. Using potfile load, the
Slater potential is not calculated but read from slater.pot (the correction
potential is instead recalculated). For spin unrestricted calculations the
corresponding files are slaterA.pot, slaterB.pot, corrctA.pot and
correctB.pot.
allows the user to specify which occupied orbital will not be included in the
calculation of correction potential: by default the highest occupied orbital is
selected. This option is useful for those systems where the HOMO of the
starting orbitals (e.g. EHT, HF) is different from the final LHF HOMO. homob
is for the beta spin.
a correlation functional can be added to the LHF potential: use func=lyp for
LYP, or func=vwn for VWN5 correlation.
For expert users
Options for the conjugategradient algorithm for the computation of the correction
potential: rmsconvergence (conjgrad conv=1.d7), maximum number of iteration
(maxit=20), output level output=03, asymptotic continuation in each iteration
(cgasy=1).
With slaterdtresh= 1.d9 (default) the calculations of the numerical integrals for the Slater potential is performed only if it changes more than 1.d9.
Asymptotic regions specification:
corrctregion R_{F }Δ_{F }
0…R_{F }  Δ_{F } : basisset correction potential
R_{F }  Δ_{F }…R_{F } + Δ_{F } : smooth region
R_{F } + Δ_{F }… + ∞ : asymptotic correction
Defaults: R_{F } = 10, Δ_{F } = 0.5
slaterregion R_{N}Δ_{N}R′_{F }Δ′_{F }
0…R_{N}  Δ_{N} : basisset Slater potential
R_{N}  Δ_{N}…R_{N} + Δ_{N} : smoothing region
R_{N} + Δ_{N}…R′_{F }  Δ′_{F } : numerical Slater
R′_{F }  Δ′_{F }…R′_{F } + Δ′_{F } : smoothing region
R′_{F } + Δ′_{F }… + ∞ : asymptotic Slater
Note: R′_{F }  Δ′_{F } ≤ R_{F }  Δ_{F }
Defaults: R_{N} = 7, Δ_{N} = 0.5, R′_{F } = 10, Δ′_{F } = 0.5
Use correctbregion and slaterbregion for the beta spin.
Selfconsistent twocomponent calculations (e.g. for spinorbit interactions) can be carried out using the module ridft . The following keywords are valid:
enforces twocomponentSCF calculations; this option is combinable with $rij,
$rik and $dft.
switches on Kramersrestricted formalism
switches on collinear twocomponent formalism (not rotational invariant)
enforces DIIS for complex Fock operator.
Relativistic allelectron calculations can be done employing the X2C, the BSS or the DKH Hamiltonian. Implemented for modules dscf and ridft.
switches on X2C calculation.
switches on BSS calculation.
switches on DKH calculation of order integer.
selects parameterization of the DKH Hamiltonian. Valid values are 1 (=default), 2,
3, 4, and 5.
Note in particular that the parametrization does not affect the Hamiltonian up to fourth order. Therefore, as long as you run calculations with DKH Hamiltonians below 5th order you may use any symbol for the parametrization as they would all yield the same results. Higherorder DKH Hamiltonians depend slightly on the chosen paramterization of the unitary transformations applied in order to decouple the Dirac Hamiltonian, but this effect can be neglected. For details on the different parametrizations of the unitary transformations see [194].
switches on local DLU approach.
selects parameterization of the local approximation. Valid values are 0 or 1. For
details on the different parametrizations see [75].
All of these keywords are combinable with $soghf.
The Periodic Electrostatic Embedded Cluster Method (PEECM) functionality provides electronic embedding of a finite, quantum mechanical cluster in a periodic, infinite array of point charges. It is implemented within HF and DFT energy and gradient TURBOMOLE modules: dscf, grad, ridft, rdgrad, and escf. Unlike embedding within a finite set of point charges the PEEC method always yields the correct electrostatic (Madelung) potential independent of the electrostatic moments of the point charges field. It is also significantly faster than the traditional finite point charges embedding.
The basic PEECM settings are defined in the $embed block. It can be redirected to an external file using $embed file=<file_name>.
Following keywords are used for the PEECM calculation setup:
Specifies the number of periodic directions. Allowed values for number are 3 for a bulk threedimensional system, 2 for a twodimensional surface slab, and 1 for a onedimensional system. Default value is 3.
Unit cell parameters in a form of six real values a, b, c, α, β, γ, where a, b, c are lengths of the appropriate cell vectors, α is the angle between vectors b and c, β is the angle between vectors a and c, and γ is the angle between vectors a and b. Default are atomic units and degrees. You can specify unit cell parameters in Å and degrees using cell ang.
Content of the unit cell, where label is the label of the point charge Content of the unit cell, where label is the label of the point charge type and x y z are corresponding Cartesian or fractional crystal coordinates. Defaults are Cartesian coordinates and atomic units. You can specify Cartesian coordinates in Å using content ang and fractional coordinates using content frac. Note that Cartesian coordinates assume that the cell vector a is aligned along the x axis, and the vector b on the xy plane.
Atomic coordinates of the piece of the crystal to be replaced by the QM cluster and surrounding isolation shell (ECPs and explicit point charges), where label is the point charge label and x y z are corresponding Cartesian or fractional crystal coordinates. Defaults are Cartesian coordinates and atomic units. You can specify Cartesian coordinates in Å using cluster ang and fractional coordinates using cluster frac.
Values of point charges (for each atomtype) , where label is the point charge label and charge specifies charge in atomic units.
Values of point charges (for each atom), where label is the point charge label and charge specifies charge in atomic units.
Note that charges and ch_list are mutually exclusive. An integer number n can also be appended to charges or ch_list to set the tolerance for charge neutrality violation to 10^{n} (default n = 5).
Additionally, the following keywords control the accuracy of PEECM calculation:
Maximum order of the multipole expansions in periodic fast multipole method
(PFMM). Default value is 25.
Electrostatic potential at the lattice points resulting from periodic point
charges field will be output if this keyword is present. Default is not to output.
Wellseparateness criterion for PFMM. Default is 3.0.
Minimum accuracy for lattice sums in PFMM. Default is 1.0d8.
The Conductorlike Screening Model (COSMO) is a continuum solvation model, where the
solute molecule forms a cavity within the dielectric continuum of permittivity epsilon that
represents the solvent. A brief description of the method is given in chapter 19. The model
is currently implemented for SCF energy and gradient calculations (dscf/ridft and
grad/rdgrad), MP2 energy calculations (RIMP2 and mpgrad) and MP2 gradients
(RIMP2), and response calculations with escf. The ricc2 implementation is described in
section 19.
For simple HF or DFT single point calculations or optimizations with standard settings,
we recommend to add the $cosmo keyword to the control file and to skip the rest of this
section.
Please note: due to improvements in the A matrix and cavity setup the COSMO energies
and gradients may differ from older versions (5.7 and older). The use_old_amat option
can be used to calculate energies (not gradients) using the old cavity algorithm of
TURBOMOLE 5.7.
The basic COSMO settings are defined in the $cosmo and the $cosmo_atoms
block.
Example with default values:
defines a finite permittivity used for scaling of the screening charges.
skips the COSMO segment statistics run and allocates memory for the given
number of segments.
skips the outlying charge correction.
All other parameters affect the generation of the surface and the construction of the A matrix:
number of basis grid points per atom
(allowed values: i = 10 × 3^{k} × 4^{l} + 2 = 12,32,42,92...)
number of segments per atom
(allowed values: i = 10 × 3^{k} × 4^{l} + 2 = 12,32,42,92...)
distance threshold for A matrix elements (Ångstrom)
distance to outer solvent sphere for cavity construction (Ångstrom)
factor for outer cavity construction in the outlying charge correction
pave intersection seams with segments
leave untidy seams between atoms
amplitude of the cavity desymmetrization
phase of the cavity desymmetrization
refractive index used for the calculation of vertical excitations and num.
frequencies (the default 1.3 will be used if not set explicitly)
in case of disjunct cavities only the largest contiguous cavity will be used
and the smaller one(s) neglected. This makes sense if an unwanted inner
cavity has been constructed e.g. in the case of fullerenes. Default is to use
all cavities.
If the $cosmo keyword is given without further specifications the default parameter are used (recommended). For the generation of the cavity, COSMO also requires the definition of atomic radii. User defined values can be provided in Ångstrom units in the data group $cosmo_atoms, e.g. for a water molecule:
If this section is missing in the control file, the default values defined in the radii.cosmo file (located in $TURBODIR/parameter) are used. A user defined value supersedes this defaults. $cosmo and $cosmo_atoms can be set interactively with the COSMO input program cosmoprep after the usual generation of the TURBOMOLE input.
The COSMO energies and total charges are listed in the result section. E.g.:
The dielectric energy of the system is already included in the total energy. OC corr denotes the outlying charge correction. The last energy entry gives the total outlying charge corrected energy in the old definition used in TURBOMOLE 5.7 and older versions. The COSMO result file, which contains the segment information, energies, and settings, can be set using: $cosmo_out file= filename.cosmo
Isodensity Cavity: This option can be used in HF/DFT single point calculations only. The $cosmo_isodens section defines the settings for the density based cavity setup (see also chapter 19). If the $cosmo_isodens keyword is given without suboptions, a scaled iosodensity cavity with default settings will be created. Possible options are:
activates the density based cavity setup. The default values of nspa and
nsph are changed to 162 and 92, respectively. This values are superseded
by the user defined nspa value of the $cosmo section. By default the scaled
density method is used. The atom type dependent density values are read
from the radii.cosmo file (located in $TURBODIR/parameter).
spacing of the marching tetrahedron grid in Å (default: 0.3Å).
use one isodensity value for all atom types (value in a.u.)
The outlying charge correction will be performed with a radii based outer cavity. Therfore, and for the smoothing of the density changes in the intersection areas of the scaled density method, radii are needed as for the standard COSMO cavity. Please note: The isodensity cavity will be constructed only once at the beginning of the SCF calculation. The density constructed from the initial mos will be used (file mos or alpha/beta in case of unrestricted calculations). Because the mos of an initial guess do not provide a good density for the cavity construction, it is necessary to provide mos of a converged SCF calculation (e.g. a COSMO calculation with standard cavity). We recommend the following three steps: perform a standard COSMO calculation, add the isodensity options afterwards, and start the calculation a second time.
Radii based Isosurface Cavity: The $cosmo_isorad section defines the radii defined isosurface cavity construction. The option uses the algorithm of the isodensity cavity construction but the objective function used depends on the cosmo radii instead of the electron density. The default values of nspa and nsph are changed to 162 and 92, respectively. This values are superseded by the user defined nspa value of the $cosmo section. The resulting surface exhibits smoother intersection seams and the segment areas are less diverse than the ones of the standard radii bases cavity construction. Because gradients are not implemented, the radii based isosurface cavity can be used in single point calculations only.
spacing of the marching tetrahedron grid in Å (default: 0.3Å).
COSMO in MP2 Calculations: The iterative COSMO PTED scheme (see chapter 19) can be used with the mp2cosmo script. Options are explained in the help message (mp2cosmo h). Both MP2 modules RIMP2 and mpgrad can be utilized. The control file can be prepared by a normal COSMO SCF input followed by a RIMP2 or mpgrad input. The PTE gradients can be switched on by using the
keyword (RIMP2 only). Again a normal SCF COSMO input followed by a RIMP2 input has to be generated. The $cosmo_correlated keyword forces dscf to keep the COSMO information needed for the following MP2 calculation and toggles on the COSMO gradient contribution.
COSMO in Numerical Frequency Calculations: NumForce can handle two types of COSMO frequency calculations. The first uses the normal relaxed COSMO energy and gradient. It can be performed with a standard dscf or ridft COSMO input without further settings. This is the right method to calculate a Hessian for optimizations. The second type, which uses the approach described in chapter 19, is implemented for ridft only. The input is the same as in the first case but Numforce has to be called with the cosmo option. If no solvent refractive index refind=REAL is given in the $cosmo section of the control file the program uses the default (1.3).
COSMO in vertical excitations and polarizabilities: COSMO is implemented in escf and will be switched on automatically by the COSMO keywords of the underlying SCF calculation. The refractive index, used for the fast term screening of the vertical excitations, needs to be defined in the cosmo section of control file (refind=REAL).
DCCOSMORS: The DCOSMORS model (see chapter 19) has been implemented for restricted and unrestricted DFT and HF energy calculations and gradients (programs: dscf/ridft and grad/rdgrad). In addition to the COSMO settings defined at the beginning of this section, the $dcosmo_rs keyword has to be set.
activates the DCOSMORS method. The file defined in this option contains
the DCOSMORS σpotential and related data (examples can be found in
the defaut potentials in the $TURBODIR/parameter directory).
If the potential file cannot be found in the local directory of the calculation, it will be searched in the $TURBODIR/parameter directory. The following σpotential files for pure solvents at 25^{∘}C are implemented in the current TURBOMOLE distribution (see parameter subdirectory):
h2o_25.pot
ethanol_25.pot
methanol_25.pot
thf_25.pot
propanone_25.pot
chcl3_25.pot
ccl4_25.pot
acetonitrile_25.pot
nitromethane_25.pot
dimethylsulfoxide_25.pot
diethylether_25.pot
hexane_25.pot
cyclohexane_25.pot
benzene_25.pot
toluene_25.pot
aniline_25.pot
The DCOSMORS energies and total charges are listed in the COSMO section of the output:
The outlying charge correction cannot be defined straight forward like in the normal COSMO model. Therefore, the output shows two corrections that can be added to the Total energy. The first one is the correction on the COSMO level (COSMO) and the second is the difference of the DCOSMORS dielectric energy calculated form the corrected and the uncorrected COSMO charges, respectively (DCOSMORS). The charges are corrected on the COSMO level only. The Total energy includes the E_{diel,RS} defined in section 19. Additionally the combinatorial contribution at infinute dilution of the COSMORS model is given in the output. The use of this energy makes sense if the molecule under consideration is different than the used solvent or not component of the solvent mixture, respectively. To be consistent one should only compare energies containing the same contributions, i.e. same outlying charge correction and with or without combinatorial contribution. Please note: the COSMORS contribution of the DCOSMORS energy depends on the reference state and the COSMORS parameterization (used in the calculation of the chosen COSMORS potential). Therefore, the DCOSMORS energies should not be used in a comparision with the gas phase energy, i.e. the calculation of solvation energies.
Many of the dscf and ridft keywords are also used by grad and rdgrad.
This keyword and corresponding options are required in gradient calculations only in special circumstances. Just $drvopt is fine, no options needed to compute derivatives of the energy with respect to nuclear coordinates within the method specified: SCF, DFT, RIDFT.
If running a DFT gradient calculation, it is possible to include the derivatives of the quadrature weights, to get more accurate results. In normal cases however those effects are marginal. An exception is numerical calculation of frequencies by Numforce, where it is strongly recommended to use the weight derivatives option. The biggest deviations from the uncorrected results are to be expected if doing gradient calculations for elements heavier than Kr using all electron basis sets and very small grids. To use the weight derivatives option, add
weight derivatives
The option
point charges
in $drvopt switches on the evaluation of derivatives with respect to coordinates of point charges. The gradients are written to the file $point_charge_gradients old gradients will be overwritten.
This module calculates analytically harmonic vibrational frequencies within the HF or (RI)DFTmethods for closedshell and spinunrestricted openshellsystems. Broken occupation numbers would lead to results without any physical meaning. Note, that RI is only used partially, which means that the resulting Hessian is only a (very good) approximation to exact second derivatives of the RIDFTenergy expression. Apart from a standard force constant calculation which predicts all (allowed and forbidden) vibrational transitions, it is also possible to specify certain irreps for which the calculation has to be done exclusively or to select only a small number of lowest eigenvalues (and eigenvectors) that are generated at reduced computational cost.
is the keyword for nondefault options of gradient and second derivative calculations.
Possibilities in case of the module aoforce are:
to read a complete Hessian from the input file $hessian and perform
only the frequency analysis
to perform an analysis of normal modes in terms of internal coordinates.
Details about this option and the effect of the printlevel (default is 0)
are given in Section 13. The effect of the keyword only is the same as
described above.
fixes the RAM memory to be used by the run (here 50MB), about 70% of available
memory should be fine, because $maxcor specifies only the memory used to store
derivatives of density and Fock matrices as well as the CPHFRHS. Default is
200MB.
sets the convergence criterion for the CPHFequations to a residual norm of 1.0d7.
Normally the default value of 1.0d5 already provides an accuracy of vibrational
frequencies of 0.01cm^{1} with respect to the values obtained for the convergence
limit.
fixes the maximum number of Davidson iterations for the solution of the
CPHFequations to a value of ten. Normal calculations should not need more than
eight iterations, but as a precaution the default value is 25.
forces the program in case of molecules with C_{1} symmetry not to use 3N  6(5)
symmetry adapted but all 3N cartesian nuclear displacement vectors. This option
may lead to a moderate speedup for molecules notedly larger than 1000 basis
functions and 100 atoms.
forces the program not to project out translations and rotations when forming a
basis of symmetry adapted molecular displacements. This option may be needed if a
Hessian is required, that contains translation and rotationcontributions, e.g. for
coupling the system with low cost methods. Output of the unprojected
hessian is done on $nprhessian; format is the same as for conventional
$hessian. Output of the corresponding eigenvalues and eigenvectors is done
analogously on $nprvibrational spectrum and $nprvibrational normal
modes.
causes the program to diagonalize a not mass weighted hessian. Output is on
$nprhessian, $nprvibrational spectrum and $nprvibrational normal modes,
because projection of rotations is not possible in this case.
This keyword allows to trace back the effects of isotopic substitution on vibrational
frequencies. The atom(s) for which isotopic substitution is to be investigated are
specified in subsequent lines of the form (atom index) (mass in special isotope),
e.g.
Sets the number of points for interpolation between the two isotopes compared by
the $isosub option to six. Default value is 21.
Keywords for the treatment of only selected nuclear displacement vectors:
CPHFiteration is done only for distortions, that are IR active.
CPHFiteration is done only for distortions, that are Raman active.
This causes a lowest Hessian eigenvalue search to be performed instead of
a complete force constant calculation. The lowest eigenvalue search consists
of the calculation of a guessHessian and macroiterations to find the
solution vector(s) for the lowest eigenvalue(s). In each macroiteration the
CPHFequations are solved for the present search vector(s). $les all 1 means
that one lowest eigenvalue for each irrep will be determined, other numbers of
lowest eigenvalues per irrep are admissible too.
Different numbers of lowest eigenvalues for different irreps are requested by e.g.
$les
a1 3
a2 all
b2 1
The convergence criterion of the Davidson iterations for the solution of the CPHFequations as well as the maximal residual norm for the lowest Hessian eigenvalue in the macroiteration are specified by $forceconv as explained above.
The maximum number of macroiterations is specified by $lesiterlimit x
with the default x=25. The maximum number of iterations for each solution of the CPHFequations is again determined by $forceiterlimit as shown above.
The convergence of the macroiterations is strongly influenced by the size
of the starting searchsubspace. Generally all guessHessian eigenvectors
corresponding to imaginary frequencies and at least two real ones are used as
starting searchsubspace. However it proved to be necessary to use even more
vectors in the case of guessHessians with very large conditioning numbers.
$hesscond 8.0d5
means that all eigenvalues with the quotient (eigenvalue)∕(max. eigenvalue)
lower than 0.00008 are added to the starting searchsubspace. Default is
1.0d4.
Triggers the generation of input files for hotFCHT (program to calculate
FranckCondonfactors by R. Berger and coworkers). See 13.4.
Save the derivative of the density matrix for subsequent use in the module
evib. See 14
Force constant calculations on the DFT level prove to be numerically reliable only with
large integration grids or if one includes the effects of quadrature weights. This is done by
default—to prevent this, insert
no weight derivatives
in $dft.
can be used to generate text output of the matrix elements of the derivative
of the Fockoperator. For bigger systems this can however generate very large
output files. See 14
to perform an escf calculation converged molecular orbitals from a HF, DFT or RIDFT calculation are needed. The HF, DFT or RIDFT method is chosen according to the $dft or $ridft keywords, specified above. It is recommended to use wellconverged orbitals, specifying $scfconv 7 and $denconv 1d7 for the groundstate calculation. The input for an escf calculation can be conveniently generated using the ex menu in define, see Section 4.
In an escf run one of the following properties can be calculated: (please note the ’or’ in the text, do only one thing at a time.)
1. RPA and timedependent DFT singlet or triplet or spinunrestricted excitation energies (HF+RI(DFT))
or
or
2. TDA (for HF: CI singles) singlet or triplet or spinunrestricted or spinflip excitation energies (HF+RI(DFT))
or
or
or
3. Twocomponent TDDFT excitation energies of Kramersrestricted closedshell systems
4. Eigenvalues of singlet or triplet or nonreal stability matrices (HF+RI(DFT), RHF)
or
or
5. Static polarizability and rotatory dispersion tensors (HF+(RI)DFT, RHF+UHF)
6. Dynamic polarizability and rotatory dispersion tensors (HF+(RI)DFT, RHF+ UHF)
$scfinstab dynpol unit
list of frequencies
where unit can be eV, nm, rcm; default is a.u. (Hartree). For example, to calculate dynamic polarizabilities at 590 nm and 400 i nm (i is the imaginary unit):
The number and symmetry labels of the excited states to be calculated is controlled by the data group $soes. Example:
will yield the 17 lowest excitations in IRREP b1g, the 23 lowest excitations in IRREP eu, and all excitations in IRREP t2g. Specify $soes alln; to calculate the n first excitations in all IRREPS. If n is not specified, all excitations in all IRREPS will be obtained.
During an escf run, a systemindependent formatted logfile will be constructed for each IRREP. It can be reused in subsequent calculations (restart or extension of eigenspace or of $rpaconv). An escf run can be interrupted by typing “touch stop” in the working directory.
The maximum amount of core memory to be allocated for the storage of trial
vectors is restricted to n MB. If the memory needed exceeds the threshold
given by $rpacor, a multiple pass algorithm will be used. However, especially
for large cases, this will increase computation time significantly. The default is
200MB.
The calculated excitation energies and corresponding oscillator strengths are
appended to a file named ’spectrum’. Possible values of unit are eV, nm and
cm^{1} or rcm. If no unit is specified, excitation energies are given in a.u.
The calculated excitation energies and corresponding rotatory strengths are
appended to a file named ’cdspectrum’. unit can have the same values as in
$spectrum.
Flag for generation of UHF start MOs in a triplet instability calculation. The
option will become effective only if there are triplet instabilities in the totally
symmetric IRREP. The optional real number e specifies the approximate
second order energy change in a.u. (default: 0.1).
Enables calculation of dipole polarizability/rotatory dispersion in the velocity
gauge. Active only for pure DFT (no HF exchange).
Enable calculation of oscillator and rotatory strength sum rules at frequencies
specified by list of frequencies in unit unit (see $scfinstab dynpol). Note that
the sums will be taken only over the states specified in $soes.
the vectors are considered as converged if the Euclidean residual norm is less
than 10^{n}. Larger values of n lead to higher accuracy. The default is a residual
norm less than 10^{5}.
Sets the upper limit for the number of Davidson Iterations to n. Default is
n = 25.
The main keyword that switches on a GW calculation. Provided that the response function is calculated setting this keyword will perform a standard G_{0}W_{0} calulation with default values for the other flags.
There are several options which can be added to the $gw keyword, the syntax is:
With the optional entries:
Default: 1. Number of orbitals to calculate gw for. It is set to the number of occupied orbitals + 5 if set smaller than the number of occupied orbitals.
Default: 0.001. Infinitesimal complex energy shift. Negative value switches to calculating at that value but extrapolating to 0 in linear approximation.
Default: qpe.dat. Output filename for the quasiparticle energies.
Default: false (not set). If added as option pure rpa responce function is calculated. If not added, the TDDFT response function is calculated and used to screen the coulomb interaction.
The keyword $rirpa allows to specify the following options,
Number of frequency integration points, n (default is 60).
HF energy calculation is skipped, (HXX = Hartree + eXact (Fock) eXchange).
Generates profiling output.
Switches on the gradients calculation for RIRPA.
Computes gradients in the direct RIMP2 limit.
Manual setting of the number of integral blocks, n, in subroutine rirhs.f. This
is for developers; the default is determined with the $maxcor.
egrad uses the same general keywords as escf and grad, see Sections 20.2.9 and 20.2.12.
The state to be optimized is by default the highest excited state specified in $soes. Note that only one IRREP can be treated at the same time in contrast to escf calculations. When the desired excited state is nearly degenerate with another state of the same symmetry it may be necessary to include higher states in the initial calculation of the excitation energy and vector in order to avoid root flipping. This is accomplished by means of the additional keyword
which explicitly enforces that nth excited state is optimized. n must not be larger than the number of states specified in $soes.
flag to compute Cartesian nonadiabatic coupling vectors between the excited state of interest and the ground state [195]. This option requires the use of weight derivatives in section dft. It is only implemented for C_{1} symmetry.
If an MP2 run is to be performed after the SCF run, the SCF run has to be done with at least
1) density convergence  $denconv 1.d7 
2) energy convergence  $scfconv 6 
The data group $maxcor adjusts the maximum size of core memory (n in MB)
which will be allocated during the MP2 run. Recommendation: 3/4 of the
actual main memory at most. If $maxcor is not found, its value is set to 200
MB.
Calculation of MP2 gradient is omitted, only MP2 energy is calculated. In
connection with this keyword you may also activate the spincomponentscaled
(SCS) MP2 proposed by Grimme
with the default values of 6/5 for pS and 1/3 for pT, which may be modified
this way:
The data group $freeze specifies frozen orbitals, in the above syntax by irreducible representations. The symmetryindependent and for standardapplications recommended syntax is
This will freeze the 5 lowest occupied and 2 highest virtual orbitals (alpha and beta count as one in UHF cases). Note that degenerate orbitals count twice (e representations), thrice (t representations) etc. In case of mpgrad frozen orbitals have to be specified manually, for rimp2 the preparation tool rimp2prep may be used to specify frozen core orbitals, frozen virtuals have to be specified manually. Note: In case of gradient calculations frozen core orbitals are regarded only by rimp2, but not by mpgrad, moreover freezing of virtual orbitals is generally not supported by mpgrad.
All essential data groups for mpgrad may be generated by the preparation tool mp2prep, apart from $maxcor (see above) these are the following:
specifies the number of loops (or ’passes’) over occupied orbitals, n, performed
in the mpgrad run: the more passes the smaller file space requirements—but
CPU time will go up.
The data group $mointunit specifies:
statistics run (estimation of disc space needed) for the adjustment of the file
sizes will be performed.
calculation of MP2 pair correlation energies.
Apart from keywords $maxcor, $mp2energy and $freeze (see above) rimp2 also needs
cross reference for the file specifying the auxiliary basis as referenced in $atoms.
We strongly recommend using auxbasis sets optimized for the corresponding
MO basis sets.
Reasonable settings for these keywords may be generated by the tool Rimp2prep. Moreover you may specify by hand:
specification of directory for scratch files; by default files are written to the
working directory; works also with capital letters (for consistency with ricc2).
avoids symmetry gymnastics in case of C_{1}symmetry, rather for debugging
enforces calculation of
, 
Enforces plotting of five largest tamplitudes as well as five largest norms of
tamplitudes for fixed pair of occupied orbitals ij. By additional integer this number
may be changed.
Enforces plotting of all eigenvalues of the MP2 density matrix.
Note that beside the keywords listed below the outcome of the ricc2 program also depends on the settings of most thresholds that influence the integral screening (e.g. $denconv, $scfconv, $scftol) and for the solution of Z vector equation with 4index integrals (for relaxed properties and gradients) on the settings for integrals storage in semidirect SCF runs (i.e. $thime, $thize, $scfintunit). For the explanation of these keywords see Section 20.2.6.
Auxiliary basis set for RI approximation. For details Section 20.2.16.
Freeze orbitals in the calculation of correlation and excitation energies. For
details see Section 20.2.16.
Print level. The default value is 1.
Specify a directory for large intermediate files (typically threeindex coulomb
integrals and similar intermediates), which is different from the directory where
the ricc2 program is started.
The data group $maxcor adjusts the maximum size of core memory in MB
which will be allocated during the RICC2 run. This keyword can be set in
define or with the Rimp2prep tool, the default is 20MB.
$maxcor has a large influence on computation times for RICC2 runs! It is
recommended to set $maxcor to ca. 75–85% of the available (physical) core
memory.
The calculated excitation energies and corresponding oscillator strengths are
appended to a file named ’spectrum’. Possible values of unit are eV, nm and
cm^{1} or rcm. If no unit is specified, excitation energies are given in a.u.
The calculated excitation energies and corresponding rotatory strengths are
appended to a file named ’cdspectrum’. unit can have the same values as in
$spectrum.
The purpose of this data group is twofold: It activates the Laplacetransformed implementation of SOSMP2 in the ricc2 module (if the sos option has been specified in $ricc2) and it provides the options to specify the technical details for the numerical Laplacetransformation.
Threshold for the numerical integration used for the Laplace
transformation of orbital energy denominators. The grid points for the
numerical integration are determined such that is the remaining root
mean squared error (RMSE) of the Laplace transformation is < 10^{conv}.
By default the threshold is set to the value of conv given in $ricc2 (see
below).
specifies the ab initio models (methods) for ground and excited states and the most important parameters and thresholds for the solution of the cluster equations, linear response equations or eigenvalue problems. If more than one model is given, the corresponding calculations are performed successively. Note: The CCS ground state energy is identical with the SCF reference energy, CCS excitation energies are identical to CIS excitation energies. The MP2 results is equivalent to the result from the rimp2 module. cis(dinf) denotes the iterative CIS(D) variant CIS(D_{∞}). The option ccsd(t) request a CCSD calculation with the perturbative triples correction, CCSD(T), and as a side result also the CCSD[T] energy will be printed.
Request the calculation of the D_{1} diagnostic in MP2 energy calculations
(ignored in MP2 gradient calculations). Note that the evaluation of the
D_{1} diagnostic increases the computational costs of the RIMP2 energy
calculation roughly by a factor of 3.
If the energy only flag is given after the cis(d) keyword, it is assumed
that only excitation energies are requested. This switches on some
shortcuts to avoid the computation of intermediates needed e.g. for the
generation of improved start vectors for CC2.
If the restart flag is set, the program will try to restart the CC2
calculations from previous solution vectors on file. If the norestart flag
is set no restart will be done. Default is to do a restart for CC2 if and only
if the file CCR0110 exists. Note: There is no restart possibility
for CCS/CIS or MP2/CIS(D).
If the hard_restart flag is set, the program will try to reuse integrals
and intermediates from a previous calculation. This requires that the
restart.cc file has been kept, which contains check sums and some other
informations needed. The hard_restart flag is switched on by default,
if the restart.cc file is present.
The conv parameter gives the convergence threshold for the CC2 ground state energy as 10^{conv}. The default value is taken from the data group $deneps.
The oconv parameter gives an additional threshold for the residual of
the cluster equations (vector function). If this parameter is given, the
iterations for the cluster equations are not stopped before the norm of the
residual is < 10^{oconv}. By default the threshold is set to oconv =conv1,
or 10× deneps if no input for conv is given.
If the norm of a vector is smaller than 10^{lindep}, the vector is assumed
to be zero. This threshold is also used to test if a set of vectors is linear
dependent. The default threshold is 10^{15}.
gives the maximum number of iterations for the solution of the cluster
equations, eigenvalue problems or response equations (default: 25).
is the maximum number of vectors used in the DIIS procedures for CC2
ground state or excitation energies (default: 10).
the maximum dimension of the reduced space in the solution of linear
equations (default: 100).
print level, by default set to 1 or (if given) the the value of the
$printlevel data group.
Fortran print format used to print several results (in particular
oneelectron properties and transition moments) to standard output.
specify wavefunction and electronic state for which a geometry
optimization is intended. For this model the gradient will be calculated
and the energy and gradient will be written onto the data groups $energy
and $grad. Required for geometry optimizations using the jobex script.
Note, that in the present version gradients are only available for ground
states at the MP2 and CC2 and for excited states at the CC2 level
and not for ROHF based openshell calculations. Not set by default.
The default model is CC2, the default electronic state the ground state.
To obtain gradients for the lowest excited state (of those included in
the excitation energy calculation, but else of arbitrary multiplicity and
symmetry) the short cut s1 can be used. x is treated as synonym for the
ground state.
the opposite–spin scaling factor cos and the same–spin scaling factor
css can be chosen. If scs is set without further input, the SCS parameters
cos=6/5 and css=1/3 are applied. This keyword can presently only be
used in connection with MP2.
the SOS parameters cos=1.3 and css=0.0 are applied. This keyword can
presently only be used in connection with MP2.
request the calculation of the D_{1} diagnostic for the ground state
wavefunction. Only needed for MP2 (see above for the alternative input
option mp2 d1diag). For all other correlated methods the D_{1} diagnostic
is evaluated by default (without significant extra costs).
calculates the secondorder corrections to the CCSD(T) energy from the
interferencecorrected MP2F12 (INTMP2F12) if $rir12 is switched
on. It can be combined either with the mp2 or the ccsd(t) methods.
In the latter case, the CCSD(T)INTF12 energy is printed. The
intcorr all keyword writes on the output all pair energies.
char=1, 2* or 2
The ansatz flag determines which ansatz is used to calculate the
RIMP2F12 ground state energy.
(Ansatz 2 is used if ansatz is absent.)
char=A, A’ or B
The r12model flag determines which approximation model is used to
calculate the RIMP2F12 ground state energy.
(Ansatz B is used if r12model is absent.)
char=F+K or T+V
The comaprox flag determines the method used to approximate the
commutator integrals [T,f_{12}].
(Approximation T+V is used if comaprox is absent.)
char=svd or cho
The cabs flag determines the method used to orthogonalize the orbitals
of the CABS basis. val is the threshold below which CABS orbitals are
removed from the calculation.
(svd 1.0d08 is used if cabs is absent.)
char=noinv, fixed or inv with flip or noflip
The examp flag determines which methods are used to determine
the F12 amplitudes. For inv the amplitudes are optimized using
the orbitalinvariant method. For fixed and noinv only the diagonal
amplitudes are nonzero and are either predetermined using the
coalescence conditions (fixed), or optimized (noinv—not orbital
invariant). If char=inv, the F12 energy contribution is computed using
all three methods. For openshell calculations noflip supresses the use
of spinflipped geminal functions.
(The fixed flip method is used if examp is absent.)
char=off or on
If char=off (default), the print out of the standard and F12
contributions to the pair energies is suppressed. The summary of the
RIMP2F12 correlation energies is always printed out.
char=LCG or R12
The corrfac flag determines which correlation factor is used for
the geminal basis. LCG requires the data group $lcg, which contains
the information regarding exponents and coefficients of the linear
combination of Gaussians.
char=off or on
The cabsingles flag determines whether or not the single excitations
into the CABS basis are computed.
The CABS singles are computed in any case if the CABS Fock matrix
elements are computed anyway for the F12 calculation (i.e., for ansatz
2 or r12model B or comaprox F+K).
char=hf, rohf, boys or pipek
The r12orb flag controls which orbitals are used for the F12 geminal basis
functions. With hf the (semi)canonical Hartree–Fock orbitals are used
(default). For ROHFbased UMP2 calculations rohf orbitals can be used,
which also implies that the $freeze data group options refer to ROHF
rather than semicanonical orbitals. For closedshell species, localised
orbitals can be used with either the Boys or PipekMezey method. For the
non(semi)canonical options, the r12orb noinv F12 energy correction is
evaluated using active occupied orbitals transformed to the same basis
as the orbitals in the geminal function.
defines the approximation to CCSDF12 which will be used if the MP2F12
calculation is followed by a CCSD or CCSD(T) calculation. The available
approximation and corresponding labels are
CCSD(F12)  ccsd(f12) 
CCSD(F12*)  ccsd(f12*) 
CCSD[F12]  ccsd[f12] 
CCSDF12b  ccsdf12b 
CCSD(2*)_{F12}  ccsd(2)_/f12 
CCSD(2)_{F12}  ccsd(2)_/f12 
It is recommended that these approximations are only used in combination with ansatz 2 and the SP approach (i.e. geminal coefficients fixed by the cusp conditions). For CCSDF12b calculations also the CCSDF12a energies are calculated as a byproduct. By default a CCSD(F12) calculation is carried out, but it is recommended that whenever appropriate the computationally more efficient CCSD(F12*) approximation is used.
In this data group you have to give additional input for calculations on excited states:
the irreducible representation.
spin multiplicity (1 for singlet, 3 for triplet); default: singlet, not needed
for UHF.
the number of excited states to be calculated within this irrep and for this multiplicity.
the number of roots used in preoptimization steps (default: npre = nexc).
the number of start vectors generated or read from file (default:
nstart = npre).
This flag switches on the calculation of oscillator strengths for excited
state—ground state transitions. Setting the parameter states=all is
mandatory for the calculation of transition properties in the present
version. The operators flag can be followed by a list of operators (see
below) for which the transition properties will be calculated. Default
is to compute the oscillator strengths for all components of the dipole
operator.
This flag switches on the calculation of oscillator strengths for excited
state—excited state transitions. Specifying the initial and final states
via istates=all and fstates=all is mandatory for the calculation of
transition properties in the present version. The operators flag can
be followed by a list of operators (see below) for which the transition
properties will be calculated. Default is to compute the oscillator
strengths for all components of the dipole operator.
require calculation of firstorder properties for excited states. For the
states option see spectrum option above; for details for the operators
input see below.
request calculation of the gradient for the total energy of an excited
state. If no state is specified, the gradient will be calculated for the lowest
excited state included in the calculation of excitation energies (Note
that only a single state should be specified; simultaneous calculation of
gradients for several states is in the present version not possible.).
convergence threshold for norm of residual vectors in eigenvalue problems is set to 10^{conv}. If not given, a default value is used, which is chosen as max(10^{conv},10^{oconv},10^{6}), where conv refers to the values given in the data group $ricc2.
convergence threshold used for preoptimization of CC2 eigenvectors is
set to 10^{preopt} (default: 3).
threshold (10^{thrdiis}) for residual norm below which DIIS extrapolation
is switched on in the modified Davidson algorithm for the nonlinear CC2
eigenvalue problem (default: 2).
If the flag leftopt is set the left eigenvectors are computed (default is
to compute the right eigenvectors, for test purposes only).
The bothsides flag enforces the calculation of both, the left and the
right eigenvectors (for test purposes only).
The oldnorm flag switches the program to the old normalization of the
eigenvectors and %T_{1} and %T_{2} diagnostics which were identical with
those used in the CC response code of the Dalton program.
In this data group you have to give additional input for the calculation of ground state properties and the solution of response equations:
This flag switches on the calculation of ground state firstorder properties (expectation values). The operators flag can be followed by a list of operators (see below) for which the firstorder properties will be calculated. Default is to compute the components of the dipole and the quadrupole moment. The option unrelaxed_only suppress the calculation of orbitalrelaxed firstorder properties, which require solution the CPHFlike Zvector equations. Default is the calculation of unrelaxed and orbitalrelaxed firstorder properties. The unrelaxed_only option will be ignored, if the calculation of gradients is requested (see gradient option below and geoopt in data group $ricc2).
requests the calculation of ground state secondorder properties as e.g. dipole polarizabilities. The operators flag has to be followed by a comma seperated pair of operators. (If more pairs are needed they have to be given with additional sop commands.) Default is to compute all symmetryallowed elements of the dipoledipole polarizability. With the freq flag on can specify a frequency (default is to compute static polarizabilities). The relaxed flag switched from the unrelaxed approach, which is used by default, to the orbitalrelaxed approach. Note that the orbitalrelaxed approach can not only be used in the static limit (freq=0.0d0). For further restrictions for the computation of secondorder properties check Chapter 10.5.
require calculation of geometric gradients. In difference to the geoopt
keyword in the data group $ricc2 this can be used to compute gradients
for several methods within a loop over models; but gradients and energies
will not be written to the data groups $grad and $energy as needed
for geometry optimizations. Note, that in the present version gradients
are only available for MP2 and CC2 and only for a closedshell RHF
reference.
convergence threshold for norm of residual vectors in linear response
equations is set to 10^{conv}. If not given in the $response data group, a
default value is used, which is chosen as max(10^{conv},
10^{oconv},10^{6}), where conv and oconv refer to the values given in the
data group $ricc2.
convergence threshold for the norm of the residual vector in the solution
of the Z vector equations will be set to 10^{zconv}.
use semicanonical formulation for the calculation of (transition)
oneelectron densities. Switched on by default. The semicanonical
formulation is usually computationally more efficient than the
noncanonical formulation. Exceptions are systems with many nearly
degenerate pairs of occupied orbitals, which have to be treated in a
noncanonical way anyway. (See also explanation for thrsemi below).
use noncanonical formulation for the calculation of (transition)
oneelectron densities. Default is to use the semicanonical formulation.
the threshold for the selection of nearly degenerate pairs of occupied
orbitals which (if contributing to the density) have to be treated in
a noncanonical fashion will be set to 10^{thrsemi}. If set to tight the
semicanonical algorithm will become inefficient, if the threshold is to
large the algorithm will become numerical unstable
threshold for preoptimizating the socalled Z vector (i.e. the lagrangian
multipliers for orbital coefficients) with a preceding RICPHF calculation
with the cbas auxiliary basis. The RICPHF equations will be
converged to a residual error < 10^{zpreopt}. Default is zpreopt=4. This
preoptimization can reduce significantly the computational costs for the
solution of the Z vector equations for large basis sets, in particular if they
contain diffuse basis functions. For calculations on large molecules with
small or medium sized basis sets the preoptimization becomes inefficient
compared to the large effects of integral screening for the conventional
CPHF equations and should be disabled. This option is automatically
disabled for ricc2 calculations based on foregoing RIJK HartreeFock
calculation.
disable the preoptimization of the Z vector by a preceding RICPHF
calculation with the cbas basis set. (Note that the preoptimization is
automatically deactivated if the ricc2 calculation is based on a foregoing
RIJK HartreeFock calculation.)
Common options for keywords in the data groups $ricc2, $response, and $excitations:
input of operator labels for firstorder properties, transition moments, etc. Currently
implemented operators/labels are
overlap (charge) operator: the integrals evaluated in the AO basis are ⟨μν⟩
dipole operator in length gauge: ⟨μr_{i}^{O}ν⟩ with i = x, y, z; the
index O indicates dependency on the origin (for expectation values
of charged molecules), which in the present version is fixed to
(0,0,0)
(all three components, individual components can be specified with
the labels xdiplen, ydiplen, zdiplen).
dipole operator in velocity gauge: ⟨μ∇_{i}ν⟩
(all three components, individual components can be specified with
the labels xdipvel, ydipvel, zdipvel).
quadrupole operator ⟨μr_{i}^{O}r_{j}^{O}ν⟩
(all six components, individual components can be specified with the
labels xxqudlen, xyqudlen, xzqudlen, yyqudlen, yzqudlen,
zzqudlen).
If all six components are present, the program will automatically give
the electronic second moment tensor (which involves only the electronic
contributions) M_{ij}, the isotropic second moment α = trM and the
anisotropy

Furthermore the traceless quadrupole moment

(including nuclear contributions) is given.
angular momentum ⟨μL_{i}^{O}ν⟩
(all three components, individual components can be specified with the
labels xangmom, yangmom, zangmom).
electronic force on nuclei ⟨μν⟩, where Z_{I} is the charge of the nucleus I and r^{I} is the position vector of the electron relative to the nucleus (all three components for all nuclei: the labels are xnef001, ynef001, znef001, xnef002, etc. where the number depends on the order in the coord file).
specification of states for which transition moments or firstorder properties
are to be calculated. The default is all, i.e. the calculations will be done
for all excited states for which excitation energies have been calculated.
Alternatively, one can select a subset of these listed in parentheses,
e.g. states=(ag{3} 1,35; b1u{1} 13; b2u4). This will select the triplet a_{g}
states no. 1, 3, 4, 5 and the singlet b_{1u} states no. 1, 2, 3 and the singlet (which is
default if no {} is found) b_{2u} state no. 4.
The specification of initial and final states for transition properties between excited
states is mandatory. The syntax is analog to the states option, i.e. either all or a
list of of states is required.
Calculate the doublesubstitutionbased diagnostics D_{2}.
Write MP2/CC2 natural occupation numbers and natural orbitals to a
file.
Calculate the error functional δ_{RI} for the RI approximation of (aibj) integrals
and its gradients with respect to exponents and coefficients of the auxiliary basis set as specified in the data group $cbas. The results are written to $egrad scaled by the factor given with the keyword $cgrad and can be used to optimize auxiliary basis sets for RIMP2 and RICC2 calculations (see Section 1.5).
define what kind of nonlinear parameters are to be optimized by relax and
specify some control variables for parameter update.
Available options are:
optimize molecular structures in the space of internal coordinates using
definitions of internal coordinates given in $intdef as described in
Section 4.1 ( default: on).
optimize molecular structures in redundant internal coordinates using
definitions of redundant internal coordinates given in $redundant. For
an optimization in redundant internal coordinates option internal has
to be switched on too, and option cartesian has to be switched off
(default: on).
optimize molecular structures in the space of (symmetrydistinct)
cartesian coordinates (default: off).
optimize basis set exponents (default=off).
Available suboptions are:
exponents of uncontracted basis functions will be optimized after
conversion into their logarithms (this improves the condition of
the approximate force constant matrix obtained by variable metric
methods and the behavior of the optimization procedure); scale
factors of contracted basis functions will not be affected by the
logarithm suboption
ALL basis set exponents will be optimized as scale factors
(i.e. contracted blocks and single functions will be treated in the
same way); if both suboptions (scale and logarithm) are given the
logarithms of the scale factors will be optimized
optimize a global scaling factor for all basis set exponents (default:
off).
define some variables controlling the update of coordinates.
Available options are:
maximum allowed total change for update of coordinates. The maximum
change of individual coordinate will be limited to dq_{max}∕2 and the
collective change dq will be damped by dq_{max}∕⟨dq∣dq⟩ if ⟨dq∣dq⟩ >
dq_{max}q
(default: 0.3)
calculate geometry update by inter/extrapolation of geometries of the
last two cycles (the interpolate option is always switched on by default,
but it is only active ANY time if steepest descent update has been chosen,
i.e. $forceupdate method=none; otherwise it will only be activated if
the DIIS update for the geometry is expected to fail)
provide a statistics output in each optimization cycle by displaying all
(the last integer, default setting by define is 5) subsequent coordinates,
gradient and energy values (default: on).
the presence of this keyword forces relax to provide informational output about the
usage of DIIS for the update of the molecular geometry.
special input related to the transformation of atomic coordinates between cartesian
and internal coordinate spaces (default: off).
Available options are:
maximum number of iterations for the iterative conversion procedure
internal → cartesian coordinates (default: 25).
convergence criterion for the coordinate conversion (default: 1.d10).
this switch activates special tasks: transform coordinates/gradients/ hessians
between spaces of internal/cartesian coordinates using the definitions of
internal coordinates given in $intdef:
available suboptions are:
the direction of the transformation is indicated by the direction of the arrow
this data group defines both the method for updating the approximate force
constant matrix and some control variables needed for the force constant
update.
Options for method:
no update (steepest descent)
Murtagh–Sargent update
Davidon–Fletcher–Powell update
Broyden–Fletcher–Goldfarb–Shanno update
combined (bfgs+dfp) update
Schlegel update
Ahlrichs update (macro option)
number of structures used
maximum number of geometries (= rank of the update procedure, for ahlrichs only)
minimum number of geometries needed to start
update
for an explanation see suboptions pulay given below e.g. ahlrichs numgeo=7 mingeo=3 maxgeo=4 modus=<gdg> dynamic
NOTES:
try to find an optimal linear combination of the coordinates of the numpul previous
optimization cycles as specified by modus (see below).
Available suboptions are:
char=<gg> or <gdq> or <dqdq> defines the quantity to be minimized
(dq = internal coordinate change).
fmode specifies the force constants to be used (only if char=<gdq> or
<dqdq>!)
fmode=static: use static force constants
fmode=dynamic: use updated force constants
real defines the threshold for the quantity g * dq∕g*dq which defines
the angle between gradient vector and coordinate change (default: 0.1).
If pulay is used in connection with a multidimensional BFGS update for
the hessian than the default is real=0.0. If > real the pulay
update for the geometry is expected to fail and will be ignored. For
example:
pulay numpul=4 maxpul=4 minpul=3 modus=<dqdq> static fail=0.2
options for $forceupdate
update only the diagonal force constants (update for offdiagonals will
be suppressed) (only active if method=ms, dfp, bfgs)
this allows to damp offdiagonal force constants by 1/real (compare
offreset, which discards offdiagonals completely). Only values > 1.0 will
be accepted. This option is active only within one relax run and will be
disabled automatically by relax. This is useful in difficult cases, where
the nondiagonal update has lead to too large nondiagonal elements of
the hessian.
reset offdiagonal force constants to zero. This option will be active for
the current optimization cycle only, i.e. it will be removed by relax after
having discarded offdiagonals!
optimization cycle specification of a maximum energy change allowed
(given in mHartree) which will be accepted using the actual approximate
force constant matrix from $forceapprox; if this energy change will be
exceeded, the force constants will be scaled appropriately
(The default: 0.0 means NO action)
scaling factor for the input hessian (default: 1.0).
lower bound for eigenvalues of the approximate hessian (default: 0.005);
if any eigenvalue drops below threig, it will be shifted to a reasonable
value defined by:
default: texttt0.005.
upper bound for eigenvalues of the hessian; if any eigenvalue exceeds
thrbig, it will limited to this value (default: 1000.0).
damp the variable metric update for the hessian by 1∕(1+ real) (default:
0.0).
specify initialization of the (approximate) force constant matrix.
Available options are:
this activates or deactivates initialization; if on has been set, relax will
provide an initial force constant matrix as specified by one of the possible
initialization options as described below and will store this matrix in
data group $forceapprox; after initialization relax resets $forceinit
to off!
provide a diagonal force constant matrix with:
available suboptions are:
this will lead to an assignment of diagonal elements (default: 1.0)).
this will lead to an assignment of initial force constant diagonals
depending on the coordinate type.
Provide individual defined force constant diagonals for
This does not work for basis set optimization. For the correct syntax of ‘fdiag=..’ see descriptions of $intdef, $global
read a cartesian (e.g. analytical) hessian from $hessian and use it as a
start force constant matrix; if $optimize internal has been set: use
its transform in internal coordinate space. If large molecules
are to be optimized, it may be necessary (large core memory
requirements!) to deactivate the numerical evaluation of the derivative of
the Bmatrix with respect to cartesian coordinates, which is
needed to transform H(cart) → H(int) exactly by specifying
no dbdx.
These keywords depend on the optimization task to be processed and are updated
by the corresponding program (i.g. SCF energy).
This data block contains nondefault specifications for the mmatrix diagonals. This
is of use if some cartesian atomic coordinates shall be kept fixed during
optimization.
Available options are:
atomic index followed by diagonal elements of the mmatrix for this atom
The scratch file ftmp allocated by relax can be placed anywhere in your file
systems instead of the working directory by referencing its pathname in this data
group as follows:
Definitions of internal coordinates and, optionally, values of internal
coordinates (val=..., given in a.u. or degrees) or force constant diagonal
elements (fdiag=...).
Cartesian coordinates and gradients calculated in subsequent optimization
cycles. Entries are accumulated by one of the gradient programs (grad, mpgrad,
rimp2, ricc2, egrad, etc.).
Basis set exponents scale factors and their gradients as calculated in subsequent
optimization cycles. Entries are accumulated by one of the gradient programs.
Global scale factors and gradients as calculated in subsequent optimization
cycles. Entries are accumulated by the grad or aoforce program.
Allows to augment internal SCF gradients by approximate increments obtained from
treatments (e.g. correlation or relativistic) on higher level. See the example
below.
Approximate force constant matrix (as needed for geometry optimization tasks).
The storage format may be specified by the
available options:
the default format is format=(8f10.5), but other 10digit f10.x formats
(e.g. x=4,6,..) are possible and will be used, after being manually specified
within $forceapprox. See the example below:
this data block contains the analytical cartesian force constant matrix (with
translational and rotational combinations projected out) as output by the
aoforce program and may be used to supply a high quality force constant
matrix $forceapprox for geometry optimizations (specifying $forceinit on
carthess, or $interconversion cartesian –> internal hessian).
either updated cartesian coordinates if a successful coordinate update has been
performed, or cartesian coordinates for input internal coordinates if only a
conversion from internal to cartesian coordinates has been performed.
updated basis set exponents, basis sets contraction coefficients or scaling
factors, if $optimize basis on has been specified.
updated global scaling factor for all basis set exponents, if $optimize global
on has been specified.
an approximate force constant matrix to be used in quasiNewton type
geometry optimizations; this matrix will be improved in subsequent
optimization cycles if one of the variablemetric methods ($forceupdate) has
been chosen. See 5.3.13 and 20.2.18.
a static (i.e. never updated) approximate force constant matrix to be used in
DIIStype geometry optimizations. It will be initialized by relax specifying:
$forceupdate pulay …modus=<dqdq> static.
The next data groups are output by relax (depending on the optimization subject) in order to control the convergence of optimization procedures driven by the shell script jobex.
real is the absolute value of the maximum component of the corresponding
gradient.
In order to save the effort for conversion of accumulated geometry and gradient data (as needed for the force constant update or the DIIS update of the geometry) to the optimization space, within which the geometry has to be optimized, one may specify the keyword
Then the relax program accumulates all subsequent coordinates and gradient
as used in optimization in this data group (or a referenced file). This overrides
the input of old coordinate and gradient data from data blocks $grad, $egrad,
…as accumulated by the grad program.
degrees
Only nondefault values are written in the control file except:
Following options are available:
Index of the Hessian eigenvector to follow for transition structure search
(transition vector). Eigenpairs are sorted in ascending order, i.e. with
increasing eigenvalues and start with index 1. The eigenpairs corresponding to
translations and rotations are shifted to the end. For minimization the value 0
has to be specified.
Method of hessian update. For minimization default is BFGS, for TS search
default is Powell and none is for no update.
Recompute the full Hessian every N’th step during a transition state search.
The default is zero and the Hessian is read in or computed in the first step only.
If the standard Hessian update methods fail, it can help to use this keyword.
Warning: This will make the calculation much more time demanding!
Freezing transition vector index.
diagonal hessian elements for diagonal Hessian guess (default: 0.5).
Maximum allowed value for trust radius (default: 0.3).
Minimum allowed value for trust radius (default: 1.0d4).
Initial value for trust radius (default tradius: radmax = 0.3).
threshold for energy change (default: 1.0d6).
threshold for maximal displacement element (default: 1.0d3).
threshold for maximal gradient element (default: 1.0d3).
threshold for RMS of displacement (RMS = root mean square)
(default = 5.0d4)
threshold for RMS of gradient (default: 5.0d4).
All values are in atomic units.
$properties specifies the global tasks for program moloch by virtue of the following options
a missing option or a option followed by the flag off will not be taken into account. The flag active may be omitted. For most of these options (with the only exceptions of trace and cowangriffin), there are additional data groups allowing for more detailed specifications, as explained below.
if moment is active you need
to compute the 0th, 1st, 2nd and 3rd moment at the reference point 0 0 0.
if potential is active you need
to compute the electrostatic potential (pot) and/or electrostatic field (fld) and/or electrostatic field gradient (fldgrd) and/or the zeroth order contribution to the diamagnetic shielding (shld) at reference point 0 0 0.
if localization is active you need $boys to perform a boyslocalization of orbitals
with orbital energies ≥thresholad=2 Hartrees; localize with respect to locxyz=x,
y and z and write resulting orbitals to lmofile= ’lmo’. At the most sweeps=10000
orbital rotations are performed. Nondefaults may be specified using the following
suboptions:
if population analyses is active you need
to perform a Mulliken population analysis. The options specify the output data:
print molecular orbital contributions to atomic s, p, d,…populations
print molecular orbital contributions to overlap populations
print atomic netto populations
print contributions of (irreducible) representations to atomic s,p,d,…populations
print contributions of (irreducible) representations to overlap populations
or
$loewdin
to perform a Löwdin population analysis (options are invalid here). A Löwdin
population analysis is based on decomposing D instead of DS in case of a
Mulliken PA.
or
to perform a population analysis based on occupation numbers (the options are not necessary and produce some output data concerning the modified atomic orbitals):
print MO contributions to occupation numbers of modified atomic orbitals (MAOs).
print all MAOs on standard output
print all MAOs to file mao.
This kind of population analysis basically aims at socalled shared electron numbers
(SEN) between two or more atoms. By default 2, 3 and 4center contributions to
the total density are plotted if they are larger than 0.01 electrons. Thresholds may
be individually chosen, as well as the possibility to compute SENs for molecular
orbitals: $shared electron numbers
orbitals
2center threshold = real
3center threshold = real
4center threshold = real
Results of this kind of PA depend on the choice of MAOs. By default, all MAOs with eigenvalues of the atomic density matrices larger than 0.1 will be taken into account. This is a reasonable minimal basis set for most molecules. If modified atomic orbitals shall not be selected according to this criterion, the data group $mao selection has to be specified
$mao selection threshold =real;
The default criterion for the selection of MAOs is the occupation number, for which a global threshold can be specified within the same line as the keyword $maoselection. If the global criterion or threshold is not desirable for some atoms, lines of the following syntax have to be added for each atom type of these.
atom symb list nmao=i method=meth threshold=r
The parameters in this definition have the following meaning:
atom symbol
list of all atoms for which this definition should apply. The syntax for this list is as usual in TURBOMOLE, e.g. 2,3,810,12
means number of MAOs to be included
means selection criterion for MAOs. meth can be occ (default), eig,
or man string, where occ denotes selection of MAOs by occupation
numbers, eig selection by eigenvalues and man allows manual selection.
In the latter case the string (max. 8 characters) appended to man serves
as nickname for the definition of the MAOs to be chosen. This nickname
is expected to appear as the leftmost word in a line somewhere within
data group $mao selection and is followed by the indices of the modified
atomic orbitals which are to be selected.
means the threshold to be applied for the selection criteria occ or eig
(default: 0.1).
Example:
option plot is out of fashion; to plot quantities on a grid, rather use $pointval in
connection with dscf, ridft, rimp2 or egrad, as described below. If nevertheless
plot is active you need
to obtain twodimensional plot data of mo 4a1g (the plane is specified by origin and two vectors with grid range and number of grid points) which is written to file 4a1g. Several plots may be obtained (#1, #2 etc.) at the same time. Use tool ’konto’ to visualize the plot.
Note: This is the oldfashioned way to plot MOs and densities. A new—and easier—one is to use $pointval, as described below.
if fit is active you need
Each line refers to all atoms, the line specifies a spherical layer of grid points around the atoms. The number of points and their distance from the van der Waals surface [Bohr] are given (the default is 1.0).
one line only, smoothing of the layers of grid points around the molecule: the real number is used to define isopotential surfaces on which the points of the layers have to lie.
One line per element has to be specified, it contains the name of the element and the van der Waals radius in [Bohr].
Properties of RHF, UHF and (twocomponent) GHF wave functions as well as those of SCF+MP2 densities or such from excited state DFTcalculations can be directly analyzed within the respective programs (dscf, ridft, mpgrad, rimp2 and egrad). In case of spinunrestricted calculations results are given for total densities (D^{α} + D^{β}) and spin densities (D^{α}  D^{β}). If not explicitly noted otherwise, in the following "D" is the SCF density, D(SCF), in case of dscf and ridft, the MP2corrected density, D(SCF)+D(MP2), for mpgrad and rimp2 and the entire density of the excited state in case of egrad. For modules dscf and ridft the analysis of properties may be directly started by calling dscf proper (or ridft proper). In case of mpgrad and rimp2 this is possible only, if the MP2 density has already been generated, i.e. after a complete run of mpgrad or rimp2.
Functionalities of analyses are driven by the following keywords.
leads to calculation of relativistic corrections for the SCF total density in case of
dscf and ridft, for the SCF+MP2 density in case of rimp2 and mpgrad and
for that of the calculated excited state in case of egrad. Quantities calculated
are expectation values < p^{2} >,< p^{4} > and the Darwin term (∑
1∕Z_{A}*ρ(R_{A})).
yields calculation of electrostatic moments arising from nuclear charges and
total electron densities. Also without setting this keyword moments up to
quadrupole are calculated, with respect to reference point (0,0,0). Supported
extensions:
By integer i; the maximum order of moments is specified, maximum and default is i=3 (octopole moments), real numbers x,y,z allow for the specification of one or more reference points.
drives the options for population analyses. By default a Mulliken PA in the basis of
cartesian atomic orbitals (CAOs) is performed for the total density (D^{α} + D^{β})
leading to Mulliken (brutto) charges and, in case of spinunrestricted calculations
also for the spin density (D^{α}  D^{β}) leading to Mulliken (brutto) numbers for
unpaired electrons. Besides total numbers also contributions from s, p, …functions
are listed separately.
Twocomponent wavefunctions (only module ridft and only if $soghf is set): In twocomponent calculations instead of S_{z} (S_{x},S_{y},S_{z}) is written to the output. Additionally a vectorfile named spinvec.txt is written, which includes the resulting spinvector for each atom in the system (also the direction).
The following modifications and extensions are supported, if the respective commands are written in the same line as $pop:
Additional information about p_{x},p_{y},p_{z} (and analogous for d and f functions) is displayed (lengthy output).
Contributions are plotted only if arising from atoms selected by list.
Contributions smaller than thrpl are not displayed (default: 0.01).
Mulliken atomic overlap matrix is displayed.
Mulliken netto populations (diagonal elements of Mulliken overlap matrix) are calculated.
Summed Mulliken contributions for a group of molecular orbitals defined
by numbers referring to the numbering obtained e.g. from the tool
eiger. Note that occupancy of MOs is ignored, i.e. all orbitals are
treated as occupied.
Mulliken contributions for single MOs defined by numbers (independent of
whether they are occupied or not). If this option is valid, one may additionally
set
to calculate a (simulated) density of states by broadening the
discrete energy levels with Gaussians and superimposing them. The
width of each Gaussian may be set by input (default: 0.01a.u.).
The resolution (number of points) may be chosen automatically
(default values are usually sufficient to generate a satisfactory
plot) or specified by hand. The output files (dos in case of RHF
wave functions, and dos_a+b, dos_ab, dos_alpha, dos_beta; for
UHF cases) contain energies (first column), resulting DOS for the
respective energy (second column) as well as s, p, dcontributions
for the respective energy (following columns).
Example:
leads to Mulliken PA (CAObasis) for each of the eleven MOs 2333, regarding only contributions from atoms 23 and 78 (results are written to standard output) and generation of file(s) with the respective simulated density of states.
to perform a natural population analyses [154]. The possible options (specified in the
same line) are
AO  must be provided, the CAO case is not implemented. 
tw=real  Threshold t_{w} to circumvent numerical difficulties in computing O_{w} 
(default: tw=1.d6). 
idbgl=integer  Debug level 
(default: idbgl=0). 
ab  For UHF cases: Print alpha and beta density results. 
short  Print only natural electron configuration and summary. 
Example:
leads to a natural population analysis (AObasis) with printing the results of alpha and beta densities (only the electron configuration and the summary) for the atoms 1,2 and 6.
To change the NMB set for atoms, one has to add a $nbonmbblock in the control file. Example:
leads to a NMB set for Ni of 4 s, 2 p and 1dfunctions and for O of 2 s and 1 pfunctions.
to perform a population analyses based on occupation numbers [155] yielding
"shared electron numbers (SENs)" and multicenter contributions. For this method
always the total density is used, i.e. the sum of alpha and beta densities in case of
UHF, the SCF+MP2density in case of MP2 and the GHF total density for
(twocomponent)GHF.
The results of such an analysis may depend on the choice of the number of modified atomic orbitals ("MAOs"), which can be specified by an additional line; without further specification their number is calculated by the method "mix", see below. Note: One should carefully read the information concerning MAOs given in the output before looking at the numbers for atomic charges and shared electron numbers.
to specify how MAOs are selected per atom.
Available options are:
a) for the way of sorting MAOs of each atom:
MAOs are sorted according to their eigenvalue (those with largest
EW finally are chosen). This is the default.
MAOs are sorted according to their occupation; note that the
number of all occupation is NOT the number of electrons in the
system. This option is kept rather for historical reasons.
b) for the determination of the number of MAOs:
A fixed number of MAOs is taken for each atom; usually this is
the number of shells up to the complete valence shell, e.g. 5 for
BNe, 6 for NaMg, etc. Exceptions are Elements Sc (Y, La), Ti
(Zr, Hf), V (Nb, Ta) for which not all five dshells are included,
but only 2, 3 or 4, respectively. This modification leads to better
agreement with partial charges calculated by an ESPfit.
All MAOs with an eigenvalue larger than <real> are chosen;
default is <real>=0.1. This and the following two options are not
valid in connection with occ.
Maximum of numbers calculated from fix and thr=0.1 is taken.
2:1 mixture of fix and thr=0.1. This choice gives best agreement
(statistical) with charges from an ESPfit. It is the default choice.
c) for additional information about MAOs:
Eigenvalues and occupations for each MAO are written to output.
Entire information about each MAO is written to output. Lengthy.
Further for each atom the number of MAOs and the sorting mode can be specified individually in lines below this keyword. Example:
atom 1,34 eig 7
leads to choice of the 7 MAOs with largest eigenvalue at atoms 1, 34.
enables the generation of localized molecular orbitals (LMOs) using Boys
localization. By default, all occupied orbitals are included, localized orbitals are
written (by default in the AObasis) to file(s) lmo in case of RHF and lalp and lbet
in case of UHF orbitals. Note, that LMOs usually break the molecular symmetry; so,
even for symmetric cases the AO (not the SAO) basis is used for the output. The
localized orbitals are sorted with respect to the corresponding diagonal element of
the Fock matrix in the LMO basis. In order to characterize these orbitals, dominant
contributions of (canonical) MOs are written to standard output as well as
results of a Mulliken PA for each LMO (for plotting of LMOs see option
$pointval).
The keyword allows for following options (to be written in the same line):
Include only selected MOs (e.g. valence MOs) in localization procedure
(numbering as available from Eiger).
maximum number of orbital rotations to get LMOs; default value is 10000
(sometimes not enough, in particular for highly delocalised systems).
lower threshold for displaying MO and Mulliken contributions (default:
0.1).
LMOs are written to file in the CAO basis (instead of AO).
PipekMezey localization is used.
fits point charges at the positions of nuclei to electrostatic potential arising from
electric charge distribution (also possible for twocomponent calculations, for
UHF cases also for spin density). For this purpose the ("real") electrostatic
potential is calculated at spherical shells of grid points around the atoms. By
default, BraggSlater radii, r_{BS}, are taken as shell radii, for each atom
the number of points is given by 1000 ⋅ r_{BS}^{2}, the total number of points
is the sum of points for each atom reduced by the number of points of
overlapping spheres. Nondefault shells (one or more) can be specified as
follows:
$esp_fit
shell i1 s1
shell i2 s2
Integer numbers i define the number of points for the respective shell, real numbers s constants added to radii (default corresponds to one shell with s=1.0).
A parameterization very close to that by Kollman (U.C. Singh, P.A. Kollman, J. Comput. Chem. 5(2), 129145 (1984)) may be obtained by
$esp_fit kollman
Here five shells are placed around each atom with r=1.4*r_{vdW } + k, k=0pm, 20pm, 40pm, 60pm, 80pm, and r_{vdW } are the vanderWaals radii of the atoms.
drives the calculation of spacedependent molecular quantities at 3D grids, planes,
lines or single points. Without further specifications the values of densities are
plotted on a threedimensional grid adapted to the molecular size. Data are deposed
to output files (suffix ṗlt) that can be visualized directly with the gOpenMol
program. In case of RHFdscf/ridft calculations you get the total density on file
td.plt, for UHFdscf/ridft calculations one gets both values for the total density
(D^{α} + D^{β}) on td.plt and the "spin density" (D^{α}  D^{β}) on sd.plt. For
mpgrad/rimp2 calculations one gets in the RHF case the total density
(D(SCF+MP2)) on td.plt and the MP2 contribution on mp2d.plt and in the UHF
case one obtains the total density (D^{α}(SCF + MP2) + D^{β}(SCF + MP2))
on td.plt, the "spin density" (D^{α}(SCF + MP2)  D^{β}(SCF + MP2))
on td.plt, and the respective MP2 contributions on files mp2d.plt and
mp2sd.plt. For egrad it is similar, just replace in the filenames mp2 by
e.
Integration of density (if absolute value greater than eps) within a sphere (origin x,y,z, radius r) is performed for
By default the origin is at (0,0,0), the radius is chosen large enough to include the whole 3D box and all contributions are regarded (eps=0).
Data different from total and spin densities are generated by following (combinable) settings (to be written in the same line as statement $pointval):
leads to calculation of electrostatic potential arising from electron densities, nuclei and—if present—constant electric fields and point charges. The densities used for calculation of potentials are the same as above; the respective filenames are generated from those of densities by replacement of the "d" (for density) by a "p" (for potential). By "pot eonly" only the electronic contribution to the electrostatic potential is calculated.
calculation of electric field. Note, that for 3D default output format (.plt, see below) only norm is displayed. Densities used are the same as above, filenames are generated from those of densities by replacement of "d" (for density) by "f" (for field).
calculation of amplitudes of MOs specified by numbers referring to the
numbering obtained e.g. from the tool eiger in the same format. The
respective filenames are selfexplanatory and displayed in the output.
Note, that also in MP2 and excited state calculations the HF/DFT
ground state orbitals are plotted (and not natural MP2/excited orbitals).
Twocomponent cases: The density of the spinors specified by numbers referring to the numbering obtained e.g. from the file EIGS are visualized. By setting the keyword minco also the amplitudes of the spinorparts are calculated, whose weights (the probability of finding the electron in this part) lie above the threshold.
calculation of amplitudes of LMOs (previously generated by $localize)
ordered by the corresponding diagonal element of the Fock matrix in the
LMO basis.
calculation of amplitudes of NMOs (previously generated by
$natural orbitals file=natural and
$natural orbital occupation file=natural
has to be set, if additionally to one of the above quantities also the density is to be computed.
calculation of the KohnSham exchangecorrelation potential. It is only valid for DFT calculations and it works for all exchangecorrelation functionals, including LHF. Note that for hybrid functionals, only the KohnSham part of the potential will be computed (the HF part is a nonlocaloperator and can’t be plotted). For GGA functional the full potential will be computed (local and nonlocal terms)
For line plots the output file is tx.vec . For UHF calculations the output files are tx.vec (alphaspin potential) and sx.vec (betaspin potential).
For a line plot the file has three columns: 1: total potential 2: local term ( or Slaterpotential for LHF) 3: nonlocal terms or Correction term for LHF
Output formats may be specified by e.g. fmt=xyz if written to the same line as $pointval. Supported are:
in case of scalars (density, (L)MO amplitudes, electrostatic potential) this format is: (x,y,z,f(x,y,z)). In case of vectors components of the vector and its norm are displayed. This format is valid for all types of grid (3D, plane, line, points, see below), it is the default format in case of calculation of values at single points. Output file suffix is .xyz.
only for 3D, default in this case. Data are written to binary files that can be directly read by gOpenMol. Note, that this output is restricted to scalar quantities; thus in case of vectors (Efield) only the norm is plotted. Output file suffix is .plt.
only for 3D. Data are written to ASCII files that can be imported by e.g. gOpenMol. Note, that this output is restricted to scalar quantities; thus in case of vectors (Efield) only the norm is plotted. Output file suffix is .map.
a format compatible with gOpenMol for visualization of vectors v. The format is x,y,z,v_{x},v_{y},v_{z}.
for planes and lines (default in these cases). In case of a line specified by α⋅ (see below) output is α,f(x,y,z) for scalars, for vectors components and norm are displayed. vectors. Analogously, in case of planes it is α,β,f(x,y,z). The output (file suffix .vec) may be visualized by plotting programs suited for twodimensional plots. A command file (termed gnuset) to get a contour plot by gnuplot is automatically generated.
only for 3D, writes out data in Cube format which can be imported by many external visualization programs.
For 3D grids nondefault boundarys, basis vector directions, origin and resolutions may be specified as follows:
Grid vectors (automatically normalized) now are (0,1,0),(0,0,1),(1,0,0), the grid is centered at (1,1,1), and e.g. for the first direction 200 points are distributed between 2 and 2.
Grids of lower dimensionality may be specified (in the same line as $pointval) by typing either geo=plane or geo=line or geo=point The way to use is best explained by some examples:
Values are calculated at a plane spanned by vectors (0,1,0) and (0,0,1) centered at (1,1,1).
Values are calculated at a line in direction (0,1,0) centered at (0,0,1). Output format as above.
Values are calculated at the two points (7.0,5.0,3.0) and (0.0,0.0,7.0).
Planeaveraged density can be computed by
The generated file td.vec will contain the quantity
 (20.1) 
The ab initio molecular dynamics (MD) program frog needs a command file named mdmaster. The interactive Mdprep program manages the generation of mdmaster and associated files. It is always a good idea to let Mdprep check over mdmaster before starting an MD run. Mdprep has onlinehelp for all menus.
In this implementation of ab initio MD, time is divided into steps of equal duration Δt. Every step, the energy and its gradient are calculated and these are used by the frog to work out the new coordinates for the next step along the dynamical trajectory. Both the accuracy of the trajectory and the total computation time thus depend crucially on the time step chosen in Mdprep. A bad choice of timestep will result in integration errors and cause fluctuations and drift in the total energy. As a general rule of thumb, a timestep Δt should be chosen which is no longer than one tenth of the shortest vibrational period of the system to be simulated.
Note that Mdprep will transform velocities so that the total linear and angular momentum is zero. (Actually, for the Leapfrog algorithm, initial velocities are Δt∕2 before the starting time).
The following keywords are vital for frog:
Number of MD time steps to be carried out. $nsteps is decreased by 1 every
time frog is run and JOBEX md stops when $nsteps reaches 0.
Number of atoms in system.
The file containing the current position, velocity, time and timestep, that is, the
input configuration. During an MD run the $current information is generally
kept at the end of the $log file.
The file to which the trajectory should be logged, i.e. the output: t=time
(a.u.);
atomic positions x,y,z (Bohr) and symbols at t;
timestep (au) Δt;
atomic symbols and velocities x,y,z (au) at t  (Δt∕2);
kinetic energy (H) interpolated at t, ab initio potential energy (H) calculated
at t, and pressure recorded at the barrier surface (atomic units, 1 au = 29.421
TPa) during the corresponding timestep;
ab initio potential energy gradients x,y,z (H/Bohr) at t.
This file can be manipulated with LOG2? tools after the MD run (Section 1.5).
The status of the MD run is a record of the action carried out during the
previous MD step, along with the duration of that step. The format matches
that of $md_action below.
Canonical dynamics is supported using the NoséHoover thermostat. This option can be enabled in Mdprep or by the following syntax:
Here, T specifies the temperature of the thermostat in K (500 K in the example) and t specifies the thermostat relaxation time in a.u. (100 a.u. in the example). It is advisable to choose the thermostat relaxation 210 times larger than the time step. Note that userdefined actions are presently not supported in canonical dynamics mode.
These are optional keywords:
Integer random number seed
Arbitrary title
To determine the trends in kinetic energy and total energy (average values and overall drifts) it is necessary to read the history of energy statistics over the recent MD steps. The number of MD steps recorded so far in each log file are therefore kept in the $log_history entry: this is updated by the program each step. The length of records needed for reliable statistics and the number of steps over which changes are made to kinetic energy (response) are specified in $ke_control.
$barrier specifies a virtual cavity for simulating condensed phases. The optional flag, angstroms, can be used to indicate that data will be entered in Ångstrøms rather than Bohr.
can be one of orth, elps, or none, for orthorhombic, ellipsoidal, or no
barrier (the default) respectively.
are the +x,y,z sizes of the cavity. In this case, an ellipsoid with a major
axis of 20Å along y, semimajor of 15Å on z, and minor of 10Å on x.
is the Hooke’s Law force constant in atomic units of force (H/Bohr) per
length unit. Here, it is 2.0H/Bohr/Ångstrøm, a bastard combination of
units.
is the effective limit to the restorative force of the barrier. For this system,
an atom at 5Å into the barrier will feel the same force as at 1.0Å.
denotes the temperature of the cavity walls in Kelvin. If the system
quasitemperature is below this setpoint, particles will be accelerated
on their return to the interior. Alternately, they will be retarded if
the system is too warm. A temperature of 0.0K will turn off wall
temperature control, returning molecules to the system with the same
momentum as when they encountered the barrier.
specifies and/or automatically generates atomic distance constraints. The optional
flag, angstroms, can be used to indicate that data will be entered in Ångstrøms
rather than Bohr.
is the convergence criterion for application of constraints. All distances
must be within +/ tolerance of the specified constraint. Additionally,
the RMS deviation of all constrained distances must be below 2/3 of
tolerance.
is the fraction of the total distance correction to be applied on each
constraint iteration.
commands frog to find the closest A atom to each atom X that is closer
than rmax and apply const. The first type line above examines each H
atom and looks for any O atoms within 1.2Å. The shortest distance, if
any, is then fixed at 0.9Å. Similarly, the second type line binds each
F to the closest C within 1.7Å, but since const=0.0, that distance is
fixed at the current value. The third type line attaches H atoms to the
appropriate nearby C, but at the current average HC distance multiplied
by the absolute value of const.
Explicitly specified constraints are listed by atom index and supercede autogenerated constraints. A positive third number fixes the constraint at that value, while zero fixes the constraint at the current distance, and a negative number unsets the constraint.
The output of frog contains the full list of constrained atom pairs and their current constraints in explicit format.
Userdefined instructions allow the user to tell frog to change some aspect of the MD run at some point in time t=real number. The same format is used for $md_status above. Here is an example:
In this example, starting from the time 2000.0a.u., velocities are to be scaled every step to keep average total energy constant. Then, from 2500.0a.u., gradual cooling at the default rate (annealing) is to occur until the time 3000.0a.u., when free Newtonian dynamics will resume.
Here are all the possible instructions:
These commands cause velocities to be scaled so as to keep the average kinetic energy (i.e. quasitemperature) or the average total energy approximately constant. This is only possible once enough information about run history is available to give reliable statistics. (Keywords $log_history, $ke_control).
At some time during the ab initio MD run the user can specify a new value for one of the dynamical variables. The old value is discarded. Single values are given by x=real number. Vectors must be read in frog format from file=file.
In Simulated Annealing MD, the temperature of a run is lowered so as to find minimumenergy structures. Temperature may be lowered gradually by a small factor each step (anneal; default factor 0.905 over 100 steps) or lowered rapidly by reversing all uphill motion (quench; default factor 0.8 each step). The cooling factors may be changed from the default using x=. Another option allows the quenching part of the run to be logged to a separate file. Alternatively, a standard nondynamical geometry optimization can be carried out in a subdirectory (relax).
Finally, this instruction turns off any previous action and resumes free dynamics. This is the default status of an MD run.
This keyword allows to carry out Tullytype fewest switches surface hopping (SH) [196]. This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [197]; for the current implementation see Tapavicza et al. 2011 [198]. In the current implementation the surface hopping algorithm only allows switches between the first excited singlet state and the ground state. However, total energies of higher excited states can be computed during the MD simulation. The proper functioning of SH has only been tested for the option
To carry out SH dynamics simulations, the keyword $surface_hopping has to be added to the control and mdmaster file. In addition several keywords are required in the control file:
needed to compute nonadiabatic couplings; this keyword requires the
use of weight derivatives in section dft
keyword needed to collect Cartesian nonadiabatic coupling vectors
along the trajectory
keyword needed to ensure dynamics starting in S_{1}
collects time integration of excitation energies along the trajectory
collects amplitudes of the adiabatic states along the trajectory
Special caution has to be taken to control problems related to conical intersections [199,200]. At geometries where conical intersections between the ground and excited state are present, DFT often exhibits singlet instabilities, which leads to imaginary excitation energies in linear response TDDFT; in this case the MD run is terminated. This problem can be circumvented by the use of the TammDancoff approximation (TDA) to TDDFT (see 7). In addition an optional keyword for the md_master file can be used:
enforces a switch to the ground state in case the S_{1}S_{0} energy gap
drops below <real> eV. As default a switch to S_{0} is enforced if the S_{1}
TDDFTTDA excitation energy becomes negative.
Often times if a switch is enforced due to a negative TDA excitation energy the potential energy surface is discontinuous and limited numerical precision of the nuclear forces may lead to a loss of total energy conservation. In this case the nuclear velocities are rescaled to obtain a conserved total energy.
In order to control the program execution, you can use the following keywords within the control file:
Switches on the calculation of the MP2 NMR shieldings. The required SCF
shielding step will be performed in the same run. This flag will be set by the
script mp2prep.
specifies the number of loops (or ’passes’) over occupied orbitals n when doing
an MP2 calculation: the more passes the smaller file space requirements—but
CPU time will go up. This flag will be set by the script mp2prep.
Scratch file settings for an MP2 calculation. Please refer to Section 20.2.16 for
a description of the syntax. This flag will be set by the script mp2prep.
Sets the convergence threshold for the shielding constant of succeeding CPHF
iterations. The unit is ppm and the default value is 0.01.
This selects the atom number for convergence check after each cphf iteration.
After this convergence is reached all other atoms are checked, too (default: 1).
have the save meaning as in dscf (see Section 20.2.6);
Since mpshift works ’semidirect’ it uses the same integral storage.
The scratch files allocated by mpshift can be placed anywhere in your file systems
instead of the working directory by referencing their pathnames in this data group.
All possible scratch files are listed in the following example:
stands for traloop start and traloop end. Each loop or pass in MP2 chemical shift
calculations can be done individually by providing the keywords $trast and $trand.
This can be used to do a simple parallelization of the run:
Create separate inputs for each traloop. Add
in the control files, number goes from 1 to the number of $traloops. Each calculation will create a restart file called restart.mpshift. To collect all steps and to do the remaining work, copy all restart files to one directory and rename them to restart.mpshift.number, add $trast 1 and $trand number_of_traloops to the control file and start mpshift.
On all systems the parallel input preparation is done automatically. Details for the parallel installation are given in Section 3.2.1. The following keywords are necessary for all parallel runs:
$parallel_platform architecture
Currently the following parallel platforms are supported:
for systems with very fast communication; all CPUs are used for the
linear algebra part. Synonyms for SMP are:
HP VClass, SP3SMP and HP S/XClass
for systems with fast communication like FastEthernet, the number of
CPUs that will be taken for linear algebra part depends on the size of
the matrices. Synonyms for MPP are:
SP3 and linuxcluster
for systems with slow communication, the linear algebra part will be
done on one single node. Synonyms for cluster are:
HP Cluster and every platform that is not known by TURBOMOLE
similar to SMP, but here the server task is treated differently: the MPI implementation on the SGIs would cause this task to request too much CPU time otherwise.
If you want to run mpgrad, $traloop has to be equal to or a multiple of the number of parallel workers.
For very large parallel runs it may be impossible to allocate the scratch files in the working directory. In this case the $scratch files option can be specified; an example for a dscf run is given below. The scratch directory must be accessible from all nodes.
For all programs employing density functional theory (DFT) (i.e. dscf/gradand ridft/rdgrad) $pardft can be specified:
The tasksize is the approximate number of points in one DFT task (default: 1000) and memdiv says whether the nodes are dedicated exclusively to your job (memdiv=1) or not (default: memdiv=0).
For dscf and grad runs you need a parallel statistics file which has to be generated in
advance. The filename is specified with
$2eints_shell_statistics file=DSCFparstat
or
$2eints’_shell_statistics file=GRADparstat
respectively.
The statistics files have to be generated with a single node dscf or grad run. For a dscf statistics run one uses the keywords:
and for a grad statistics run:
maxtask is the maximum number of twoelectron integral tasks,
maxdisk defines the maximum task size with respect to mass storage (MBytes) and
dynamic_fraction is the fraction of twoelectron integral tasks which will be allocated
dynamically.
For parallel grad and rdgrad runs one can also specify:
This means that the density matrix is computed by one node and distributed to the other nodes rather than computed by every slave.
In the parallel version of ridft, the first client reads in the keyword $ricore from the control file and uses the given memory for the additional RI matrices and for RIintegral storage. All other clients use the same amount of memory as the first client does, although they do not need to store any of those matrices. This leads to a better usage of the available memory per node. But in the case of a big number of auxiliary basis functions, the RI matrices may become bigger than the specified $ricore and all clients will use as much memory as those matrices would allocate even if that amount is much larger than the given memory. To omit this behavior one can use:
specifying the number of MBs that shall be used on each client.
For parallel jobex runs one has to specify all the parallel keywords needed for the different parts of the geometry optimization, i.e. those for dscf and grad, or those for ridft and rdgrad, or those for dscf and mpgrad.