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. Complete control files are provided in Chapter 21. An alphabetical list of all keywords is given in the index.
The control file and other input files references therein 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.
Most postSCF programs (aoforce, ccsdf12, egrad, escf, evib, mpgrad, mpshift, pnoccsd, ricc2, rirpa) in TURBOMOLE read the data group
to control the amount of dynamically allocated memory. The integer number (above “500”) one can (optionally) give a SI prefix (kb, mb, gb) or an IEC prefix (kib, mib, gib) to specify the unit. If no unit is specified mib (binary megabytes, i.e. 1024^{2} bytes) are used.
In addition one can specify a reference for parallel calculations:
If not given the specified memory will be used per thread. Note that in release versions before V7.2 the memory was specified per process. For sequential calculations this suboption has no effect.
If $maxcor is not given 500 MiB per thread will be used as limit for dynamically allocated memory.
Note:
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.
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.5). 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 765 and for definition
of the functionals the section 6.2 on page 250.
Example (default input):
Specification of the spherical grid (see section 20.2.7 on page 765). 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 770). For the definition of gridtype 1 – 3, please refer to Eq. (16), (17) and, (19) in Ref. [232].
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 [233],
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. [232].
In addition, the reduction of spherical grid points near nuclei is suppressed, i.e. fullshell on is set (see page 741).
Note: the keyword radsize integer overrules the setting of incr; for more information, see p. 736.
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. 750 and p. 756).
Example (Recommendation):
—for developers only—
Radial grid points have a linear scaling parameter ξ, see Eq.(16)–(19) and Table 1 in
Ref. [232]. 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. [232]. 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. 771.
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.7 on page
750).
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.12
on page 811.
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:
Within the seminumerical integration scheme employed for local hybrid functionals,
S and Pjunctions feature prescreenings of the analytical twocenter integrals with
respect to shell pair overlap and shell pair conjunction through the density matrix,
respectively. This allows to neglect certain shell pairs and thus accelerates the
calculation. sjunc and pjunc set the corresponding thresholds to 10^{integer}.
Smaller values increase the number of neglected shell pairs and thus reduce the
calculation cost, but also decrease the accuracy of the results. Note that for gradient
and scf calculations the Pjunctions also depend on the grid weights. If
larger grids are employed seeking more accurate results, the threshold for
Pjunctions should thus be adjusted (increasing the value to pjunc 6). The
defaults settings are sjunc 6 and pjunc 6 for grad, rdgrad, and ridft.
Too small values for sjunc and pjunc might result in nonconverging
calculations or deteriorated results. This holds particularly for excited state
calculations. For escf the default values are therefore sjunc 8 and pjunc
8.
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: 300 K) is reduced at each SCF cycle by the factor tmfac (default: 1.0) until it reaches the value tmend (default: 300 K). Note that the default values lead to occupation numbers calculated at a constant T = 300 K. 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. 50 K.
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.
The option addTS adds the entropic contribution of the smearing to the total energy which is important for very high temperatures only.
Finally the option noerf will use the full correct Fermi statistics rather than the term given above.
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 18).
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.
This keyword enables the use of the maximum overlap method in unrestricted ridft
calculations to access excited states selfconsistently. [234]
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.
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 PARALLEL PROCESSING
Options dscf parallel, grad, mpgrad 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 [235].
switches on local DLU approach.
selects parameterization of the local approximation. Valid values are 0 or 1. For details on
the different parametrizations see [92].
All of these keywords are combinable with $soghf.
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)
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 17.2. 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 17.2.
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. If the option
ion is added to the same input line, the scaling factor for ions f(ε) = with
x = 0 will be used. Alternatively the x value can be set by adding ion=x with x
as a real value.
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 17.2). 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 17.2) 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 17.2, 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 17.2) 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 17.2. 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.
riper shares most of the relevant keywords of the dscf and ridft modules. The $dft data group (see 20.2.7) and auxiliary basis sets defined using the keyword $jbas are always required.
For periodic calculations two additional keywords are necessary:
Specifies the number of periodic directions: n = 3 for a 3D periodic bulk solid, n = 2
for a 2D periodic surface slab and n = 1 for a 1D periodic system. The default value
is 0 for a molecular system.
Specifies the unit cell parameters. The number of cell parameters depends on the
periodicity of the system:
For 3D periodic systems six unit cell parameters a, b, c, α, β and γ need to be provided. Here, a, b and 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. riper assumes that the cell vectors a and b are aligned along the x axis and on the xy plane, respectively.
For 2D periodic systems three surface cell parameters a, b and γ have to be provided. Here, a and b are lengths of the appropriate cell vectors and γ is the angle between a and b. riper assumes that the cell vectors a and b are aligned along the x axis and on the xy plane, respectively.
For 1D periodic systems only one parameter specifying the length of the unit cell has to be provided. riper assumes that periodic direction is along the x axis.
By default, the parameters have to be provided in atomic units and degrees. Alternatively, Å can specified using $cell angs.
Alternatively, lattice vectors can be provided. The number of cell parameters depends
on the periodicity of the system:
For 3D periodic systems three (threedimensional) lattice vectors need to be provided.
For 2D periodic systems two (twodimensional) lattice vectors have to be provided. riper assumes that the lattice vectors are aligned on the xy plane.
For 1D periodic systems only one parameter specifying the length of the lattice vector has to be provided. riper assumes that periodic direction is along the x axis.
By default, lattice vectors have to be provided in atomic units. Alternatively, Å can be specified using $lattice angs.
Optionally, for periodic systems a k points mesh can be specified:
nkpoints n_{1} n_{2} n_{3}
Specifies components along each reciprocal lattice vector of a Γ centered mesh of k points.
In 3D periodic systems each k point is defined by its components k_{1}, k_{2} and k_{3} along the
reciprocal lattice vectors b_{1}, b_{2} and b_{3} as
 (20.1) 
For 2D periodic systems k_{3} = 0. In case of 1D periodicity k_{3} = 0 and k_{2} = 0. The three components k_{j} (j = 1,2,3) of k are given as
 (20.2) 
with n_{j} (j = 1,2,3) as integer numbers. For 3D periodic systems n_{1}, n_{2} and n_{3} have to be specified. The component n_{3} can be omitted for 2D periodic systems. In case of 1D periodicity only the n_{1} has to be specified.
In addition, the options kptlines and recipr within the $kpoints group can be used to obtain band structure plots:
Specifies the number of reciprocal space lines for band structure plots.
Specifies lines in the reciprocal space using three real numbers defining the start point of the line and three real numbers defining its end point. Finally, the number of k points along the line is given as an integer number.
In the following example band energies are calculated along four lines, as specified by the keyword kptlines 4. Each line definition starts in a new line with the keyword recipr, followed by three real numbers defining the start point of the line and three real numbers defining its end point. Finally, the number of k points along the line is given as an integer number. Thus, the first line starts at the point (0.500 0.500 0.500), ends at (0.000 0.000 0.000) and contains 40 k points.
The calculated band structure is written to the file bands.xyz. Each line of the file contains five real numbers: the coordinates k_{1}, k_{2} and k_{3} of the k point, its length k and the corresponding band energy ϵ_{nk}.
Calculation of energy first derivatives with respect to the lattice vectors can be requested using the keyword $optcell. The actual values and gradients of lattice vectors are written to the control file in the data group $gradlatt. Alternatively, another file name can be specified using $gradlatt file=.
Keywords within the section $riper can be used to control precision and parameters of algorithms implemented in riper. If the $riper group is absent, the following default values are used:
The following options are available:
Threshold for integrals neglect and for differential overlap when screening basis functions pairs. Probably never needs to be changed.
Flag for energy calculation only, no gradients.
If set to on charge projection of the auxiliary electron density [111] is performed for molecular systems during calculation of the Coulomb term. The charge projection constraints the charge of auxiliary density exactly to the number of electrons in the system. It is required for periodic systems, otherwise the Coulomb energy would be infinitely large. For molecular systems charge projection leads to a slight increase of the RI fitting error. It may be useful in some cases but we have so far not identified any.
Forces orthonormalization of orbital coefficients every northol SCF iteration.
If set to on full diagonalization of the Coulomb metric matrix [111] is performed and used to solve density fitting equations. When diffuse auxiliary basis functions are used the default Cholesky decomposition of the Coulomb metric matrix may fail due to small negative eigenvalues. In this case the slower method based on a full diagonalization of the metric matrix is necessary.
If pqmatdiag is used pqsingtol sets threshold for neglect of small eigenvalues of the Coulomb metric matrix.
Maximum lmoment of multipole expansions used for calculation of the Coulomb term. The default value hardly ever needs to be changed.
Target number of charge distributions per lowest level box of the octree [107]. The default value hardly ever needs to be changed.
Sets the wellseparateness criterion [107]. Octree boxes with centers separated more than sum of their lengths times 0.5×wsicl are considered as wellseparated. The default hardly ever needs to be changed. wsicl makes sense only for values ≥ 2.0. For wsicl ≤ 3.0 increasing lmaxmom may be necessary for reasonable accuracy.
Precision parameter used to determine basis function extents [107].
If set to on, an additional acceleration method employing local multipole expansions is used. For periodic systems this leads to a significant speedup of calculations, especially for small unit cells and/or diffuse basis functions. Default value is off and on for molecular and periodic systems, respectively.
For locmult set to on the order of local multipole expansions is increased by locmomadd. The default value probably never needs to be changed.
If set to on the lowmemory iterative density fitting (LMIDF) scheme is used for solving the RI equations [110] using the preconditioned conjugate gradient (PCG) method. It is implemented for molecular systems only. Default value is off.
If lpcg is used, lcfmmpcg specifies whether the CFMM is applied for evaluation of the matrixvector products needed for the PCG solver. Not employing CFMM speeds up the calculations but significantly increases memory demand since the full Coulomb metric matrix has to be stored. Default value is on.
Maximum lmoment of multipole expansions for calculations of Coulomb interactions within the PCG algorithm. It should be set to the same or larger value than lmaxmom.
Sets the threshold parameter controlling accuracy of the PCG solver (see [110] for details). Default value is 1.0 ⋅ 10^{9}. For lowerprecision calculations it can be set to 1.0 ⋅ 10^{8} but values larger than 1.0 ⋅ 10^{7} are not allowed as these lead to large errors in Coulomb energies and occasionally to SCF convergence problems.
char = at, ss or sp
Specifies the type of preconditioner used in the PCG algorithm. Three
types of preconditioners are implemented and are defined explicitly in
Sec. 7. The sp preconditioner is a default one performing consistently
the best among the preconditioners considered. The at preconditioner
is less efficient in decreasing the number of CG iterations needed
for convergence. However, it has negligible memory requirements
and more favorable scaling properties, albeit with a large prefactor.
The ss preconditioner represents a middle ground between the sp
and at preconditioners both in terms of the efficiency and memory
requirements.
Width of the Gaussian smearing in hartree. See ref. [112] for more information. Note, that [112] uses eV as the unit for the width of the Gaussian smearing.
Specifying desnue along with sigma forces occupancy leading to the number of unpaired electrons equal to desnue.
Keywords within the section $pointvalper can be used to evaluate quantities (electron density/molecular orbitals) for visualization on grid points:
The following options are available:
Specifies output format. Currently .plt, .cub, .xyz and .upt extensions are supported (see 7.3.8).
If present, total (and spin for UHF) density is calculated.
Specifies the number of plotted orbitals. For further details see subsection 7.3.8.
Number of unit cell images n1, n2 and n3 in the periodic directions a, b and c, respectively, for which plot data is generated.
Number of grid points n1, n2 and n3 along each periodic direction. If not specified, value 100 is used for each n.
Specifies the distance real in bohr around the system for which plot grid is generated in aperiodic directions. Default value is 5 bohr.
Number of grid points stored in one octree box during density calculations. Default value is 50. For very large systems or high resolution it may be necessary to increase this parameter to avoid memory allocation problems.
Only valid for plt output format. This format uses orthogonal grids. Therefore, for nonorthogonal unit cells grid data is generated for a rectangular box that contains the supercell (unit cell and its periodic images). By default, the values at grid points outside of the supercell are set to zero. For strongly nonorthogonal systems this may lead to large files. The option full switches off the zeroing of values on grid points outside the supercell.
To calculate a simulated density of states (DOS) keyword $dosper is necessary:
The following options are available:
The width of each Gaussian, default value is 0.01 a.u.
Lower/upper bounds for energy in DOS calculation.
Scaling factor for DOS (total and s, p, ... contributions).
Resolution (number of points).
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 14. The effect of the keyword only is the same as described above.
fixes the amount of core memory to be used for dynamically allocated arrays (here 50 MiB),
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 right
hande side vectors for the CPHF equations. For further details see subsection
20.2.2.
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.01 cm^{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.
Calculation of VCD rotational strenghts after preceding mpshiftrun.
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 14.5.
Save the derivative of the density matrix for subsequent use in the module evib. See
15
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 15
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/TDA excitation energies of Kramersrestricted closedshell systems
or
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 200 MB.
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} calculation 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: qpenergies.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 main keyword that switches on a BSE calculation. Provided that the response function is calculated setting this keyword will perform a standard BSE calculation with default values for the other flags. The following optional entries can be added to the $bse keyword:
Default: Not set. If set, the matrices A and W, respectively, are constructed using KS/HF orbital energies instead of GW quasiparticles energies.
Default: qpenergies.dat. This option defines the file name of GW quasiparticle energies which are read in.
Default: Not set. Compute the auxiliary matrix for the static screened interaction iteratively instead of by Cholesky decomposition.
Default: 1.0d12. Threshold for convergence in case of the iterative computation of the static screened interaction.
Default: 100. Maximum number of iterations in case of the iterative computation of the static screened 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.12 and 20.2.15.
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 [236]. 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. For further details see subsection 20.2.2.
Calculation of MP2 gradient is omitted, only MP2 energy is calculated.
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. Note: In case of gradient calculations frozen core orbitals are not regarded 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.
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.7.
Auxiliary basis set for RI approximation.
Freeze orbitals in the calculation of correlation and excitation energies. For details see
Section 20.2.18.
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 which will be
allocated during the ricc2 calculation. $maxcor has a large influence on computation
times. It is recommended to set $maxcor to ca. 75–85% of the available physical core
memory. For further details see subsection 20.2.2.
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_{∞}).
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 ground state energy for the iterative coupledcluster methods 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 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.
If no_f12metric is selected the coulomb metric is used in the density fitting
scheme to calculate the four index integrals over the operators f_{12}, f_{12}g_{12}, f_{12}^{2}
and f_{12}r_{12}. If f12metric is selected the operator’s own metric is used. The
default for the ricc2 program is no_f12metric, while the pnoccsd program
can only be used with f12metric, where it is therefore the default.
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.
request the calculation of twophoton transition moments between ground and
excited states. It recognizes the optional suboptions states for specifying
the states for which the twophoton transition moments should be computed,
operators for specifying a list of pairs of oneelectron transition operators,
and freq for a list of frequencies for one of the photons. If these suboptions
are specified, twophoton transition moments are computed per default for all
requested singlet states, dipole operators (in the length representation) and
photon energies which are equal 1/2 of the transition energies.
request the calculation of derivatives of transition moments between ground
and excited states. It recognizes the optional suboptions states for specifying
the states for which the twophoton transition moments should be computed,
operators for specifying a list of pairs of oneelectron transition operators. For
phosphorescence lifetimes operators should be set to (diplen,soc) and freq
to 0.0d0.
requests the 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. The simultaneous calculation
of gradients for several state is 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).
The ccsdf12 program uses a subset of the data groups of the ricc2 program. For F12 calculations it uses in addition the data groups $rir12, $lcg, $cabs, and $jkbas. See the respective sections for further details.
The ccsdf12 program recognizes in the $ricc2 data group the following options to specify wavefunction methods:
The options mp3, mp4, and cccsd request, respectively, MP3, MP4 and CCSD calculations. 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.
In the data group $rir12 it recognizes ccsdapprox as additional option:
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.
Note that beside the keywords listed below the outcome of the pnoccsd program also depends on the settings of most thresholds that influence the integral screening (e.g. $denconv, $scfconv, $scftol). For the explanation of these keywords see Section 20.2.7.
Auxiliary basis set for RI approximation.
Freeze orbitals in the calculation of correlation and excitation energies. For details see
Section 20.2.18.
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 program
is started.
The data group $maxcor adjusts the maximum size of core memory which will be
allocated during the pnoccsd run. $maxcor has a large influence on computation
times. It is recommended to set $maxcor to ca. 75–85% of the available physical core
memory. For further details see subsection 20.2.2.
Only needed for test purposes. It sets the accuracy for the numerical Laplacetransformation to 10^{conv}. For OSVPNOMP2 this threshold is only used for the approximate amplitudes from which the OSVs are computed when the iterative algorithm is enabled. The default value is 10^{1}. Decreasing it has little influence on the accuracy of the final results but slows down the OSV generation.
specifies the ab initio model (method). The current release version is restriced to MP2 which is also the default model.
specifies the localization method; possible choices are boys for FosterBoys, pm for PipekMezey, ibo for intrinsic bond orbitals (IBOs) and none for canonical (deprecated, only meant for testing) orbitals. Default are FosterBoys orbitals, IBOs are recommended if seperation of σ and πtype orbitals is important. (PipekMezey orbitals become very delocal with diffuse basis sets.)
switch between to algorithms for the generation of OSV coefficients. Possible choices are full for an (^{4}) scaling direct diagonalization (cmp. [237]) and davidson for an (^{3}) scaling iterative diagonalization (cmp. [20]). davidson is the default choice and recommended.
the maximal dimension of trial vectors in the iterative OSV generation (prepno davidson). For PNOMP2 the default is 800. The dimension is bounded by the number of active virtual orbitals and except for small systems a much smaller value as the number of virutals is sufficient. Smaller dimensions increase the performance, but than the iterative scheme might not converge and the program must be restarted with an adjusted dimension. For PNOMP2F12 the Default is 4000 since a larger reduced space is required to construct the OSX (cmp. [237]). Usually there is no need to touch this parameter.
specifies the PNO truncation threshold. Default value: 10^{7}.
specifies the OSV truncation threshold. If no given tolosv is set to 0.1×tolpno, which is the recommended value.
specifies the threshold for selecting orbital and pairspecific auxiliary basis sets for the local RI approximation. If not given tolri is set to 10^{7∕12} ×, which is the recommended value.
specifies the energy threshold for selecting the significant pairs. If not given tolpair is set to (0.1 ×tolpno)^{2∕3}
enables (on) or disables (off) for F12 calculations the use of OPNOs for the occupied orbital spaces in the projectors for the threeelectron integrals. Default: on
sets for F12 calculations the truncation thresholds for the complementary auxiliary PNOs for the virtual spaces (CAPNOs) in the projectors for the threeelectron integrals. If not specified, default values are calculated from the threshold tolpno.
sets for F12 calculations the truncation thresholds for the orbital specific complementary auxiliary virtuals from which the CAPNOs are generated. If not specified, default values are chosen as 0.1 × the tolcapno thresholds, which is the recommended choice.
sets for F12 calculations the truncation thresholds for the selection of PNOs for the occupied orbital spaces (OPNOs) in the projectors for the threeelectron integrals. If not specified, default values are calculated from the threshold tolpno.
set for F12 calculations the truncation thresholds for the orbital specific auxiliary occupied orbitals from which the OPNOs are generated. If not specified, default values are chosen as 0.1 × the tolopno threshold, which is the recommended choice.
sets for F12 calculations the pair specific projector ^{A}Q_{ij} for the threeelectron
integrals. Possible choices are:
projector 1 for 1  _{1} _{2}  _{1} _{2}  O_{1} _{2}^{′′} _{1}^{′′}O_{2}
projector 2 for 1  _{1} _{2}  _{1} _{2}  _{1}V _{2}^{′′} V _{1}^{′′} _{2}
projector 3 for 1  _{1} _{2}  _{1} _{2}  _{1} _{2}^{′′} _{1}^{′′} _{2}
With projector 3 ^{A}Q_{ij} becomes independend of the system size and is
therefore the default. But note that for the (^{4}) scaling algorithm this can
not be exploited and here projector 1 is more suited.
The conv parameter gives the convergence threshold for the ground state energy as 10^{conv}. The default threshold is 10^{7}.
The oconv parameter gives an additional threshold for the residual of the ground state equations as < 10^{oconv}. The default threshold is 10^{3}.
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 prePNOs is linear dependent. The default threshold is 10^{12}.
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 the ground state equations (default: 10).
not used in the current release.
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.
the SOS parameters cos=1.3 and css=0.0 are applied.
for the description of this data group see Sec. 20.2.19.
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.22.
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.01 a.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 [213]. 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 [214] 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 AO basis) 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.
FosterBoys localization is used.
Intrinsic Bond Orbitals are used.
Compute LMO centers, i.e. the diagonal matrix elements of the position
operator _{i} = ⟨ii⟩
Compute LMO spreads defined as σ_{i} = = ,
where _{i} is the LMO center
sort LMOs such that LMOs with spatially close centers are also close in
index space, by default the LMOs are sorted according to the orbital energy
expectation values
print timings for localization procedures
triggers the generations of a wfn file. It can be used in dscf/ridft singlepoint calculations or
in ricc2/egrad gradient calculations.
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}.
a format which writes out all data in extended precision as ASCII file. Useful to import in general 3D visualization tools or spread sheet programs. Selfexplaining format (open one in an editor to see its definition).
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.
In particular for 2D plots it is often useful to shift and rotate the grid plane such that some particular atoms are located in the plot plane. This can be achieved with the orient option, which accepts as additional input a list of atoms, e.g. the atoms 45 and 9 in the following example:
This will shift the origin of the plot into the center of mass of the specified set of atoms and align the grid axes with their principal axes: If two or more atoms are specified which lie on one line grid axis 1 is rotated into this line. If the specified subset consists of three atoms located in a plane the grid axis 1 and 2 are rotated into this plane. If more than three atoms are specified the grid axis 1 and 2 are rotated into a plane which fitted to the position of all specified atoms. With the rotate option the grid can be rotated around a grid axis. It accepts as input the rotation angle the index (1/2/3) of the grid axis around which the grid will be rotated.
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.3) 
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.0 H/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.0 K 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.0 a.u., velocities are to be scaled every step to keep average total energy constant. Then, from 2500.0 a.u., gradual cooling at the default rate (annealing) is to occur until the time 3000.0 a.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) [238]. This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [239]; for the current implementation see Tapavicza et al. 2011 [240]. 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 [241,242]. 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 8). 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.18 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.7);
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.
U_{ai}^{Bβ} will be written to file for a subsequent calculation of VCD rotational strengths by
aoforce.
Note: All options in this section are obsolete since the task distribution is done automatically — but they are kept for compatibility reasons.
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 can be used for parallel runs, but are all optional:
$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
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 DFT with dscf/grad$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, a parallel statistics file which has to be generated in advance can be used.
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.