20.2 Format of Keywords and Comments

TURBOMOLE input is keyword-directed. 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:


20.2.1 Keywords for System Specification

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.

$symmetry d4h

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 (C3, C3h, T, etc). Use a lower symmetry group in this case.



cu 1-4                                                             \  
   basis =cu ecp-18 arep                                           \  
   jbas  =cu ecp-18                                                \  
   ecp   =cu ecp-18 arep  
se 5-6                                                             \  
   basis =se ecp-28 arep dzp                                       \  
   jbas  =se ecp-28                                                \  
   ecp   =se arep

note the backslash \ : this is necessary. For each type of atom, one has to specify

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.

$pople char

This data group specifies the number of Cartesian components of basis functions (i.e. 5d and 7f in AO-Basis, 6d and 10f in CAO-Basis) 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).


$closed shells

Specification of MO occupation for RHF, e.g. 

 a1g     1-4                                    ( 2 )  
 a2g     1                                      ( 2 )

$open shells type=1

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:

b2g     1                                      ( 1 )  
b3g     1                                      ( 1 )  
$roothaan         1  
a = 1      b = 2


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.

$alpha shells  
 a1g     1-4                                    ( 1 )  
 a2g     1                                      ( 1 )  
$beta shells  
 a1g     1-4                                    ( 1 )  
 a2g     1                                      ( 1 )

The specification of MO occupation for UHF, $uhf overwrites closed-shell occupation specification.

20.2.2 Keyword for the General Memory Specification

Most post-SCF programs (aoforce, ccsdf12, egrad, escf, evib, mpgrad, mpshift, pnoccsd, ricc2, rirpa) in TURBOMOLE read the data group

$maxcor 500 mib per_core

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. 10242 bytes) are used.

In addition one can specify a reference for parallel calculations:

indicates that the memory is specified per process
indicates that the memory is specified per computer node, i.e. shared by all processes of the calculation that run on the machine with the same host name
indicates that the memory is specified per core, i.e. per thread
means the specified memory should be devided by all processes and nodes

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 sub-option has no effect.

If $maxcor is not given 500 MiB per thread will be used as limit for dynamically allocated memory.


20.2.3 Other General Keywords

$operating system unix  
$lock off  
$suspend off

The four keywords above are set by define, but are not necessary.

$statistics dscf
$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.

$actual step dscf

means current step. Keyword and data group (as e.g. dscf) is set by every program and removed on successful completion.

$last step relax

Keyword and data group (as e.g. relax) set by every program on successful completion.

General file cross-references:

$coord              file=coord  
$intdef             file=coord  
$user-defined bonds file=coord  
$basis         file=basis  
$ecp           file=basis  
$jbas          file=auxbasis  
$scfmo           file=mos  
$uhfmo_alpha     file=alpha  
$uhfmo_beta      file=beta  
$natural orbitals           file=natural  
$natural orbital occupation file=natural  
$energy        file=energy  
$grad          file=gradient  
$forceapprox   file=forceapprox

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.

$coord (and $intdef and $userdefined bonds)

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.


auxiliary (fitting) basis for the Coulomb terms in ridft.

$scfmo, $uhfmo_alpha, $uhfmo_beta

MO vectors of SCF or DFT calculations for RHF or UHF runs.

$natural orbitals, $natural orbital occupation

keywords and data groups set by unrestricted dscf or ridft runs. Contain natural MO vector and orbital occupation.

$energy, $grad

energies and gradients of all runs, e.g. for documentation in a geometry optimizations.


approximate force constant for geometry optimizations.

20.2.4 Keywords for redundant internal coordinates in $redund_inp

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:

iprint n

print parameter for debug output: The larger n is, the more output is printed n  0,n  5 (default: 0)

metric n

method for generating and processing of redundant internal coordinates
n   - 3,n  3,n  0 (default: 3)

Values for the metric option:

n = 1

 “Delocalized Coordinates”
The BmBt matrix is diagonalized for the complete set of redundant internal coordinates, matrix m is a unit matrix.

n = -3

Delocalized Coordinates obtained with a modified matrix m, the values of m can be defined by user input (see below).

n = -1

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

n = -2

Very simular to the n = 1 option, but for the remaining cage delocalized coordinates with modified matrix m are defined as for n = -3.

n = 2

 “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 BBt matrix is diagonalized for each block separately after the columns of B were orthogonalized against the columns of B of the the preceding blocks.

n = 3

 “Generalized natural coordinates”
Natural internal coordinates are defined first, for the remaining cage decoupled coordinates are defined.

type r

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

Types of internal coordinates for the definition of m

The matrix m is assumed to be a diagonal matrix. For each type of coordinate a different value for the force constants mii 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)

20.2.5 Keywords for Module uff

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:

              1         1          0 ! maxcycle,modus,nqeq  
         111111                      ! iterm  
       0.10D-07  0.10D-04            ! econv,gconv  
           0.00  1.10                ! qtot,dfac  
       0.10D+03  0.10D-04       0.30 ! epssteep,epssearch,dqmax  
             25      0.10       0.00 ! mxls,dhls,ahls  
           1.00      0.00       0.00 ! alpha,beta,gamma  
              F         F          F ! transform,lnumhess,lmd

The explanation of the variables are as follows:


number of max. optimization cycles (maxcycle=1: single-point 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.

econv, gconv

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 deepest-descent-step will be done.


if the norm of the gradient is smaller than epssearch, no line-search step will be done after the Newton step.


max. displacement in a.u. for a coordinate in a relax step.

mxls, dhls, ahls

parameters of linesearch:


start value




number of energy calculations

alpha, beta, gamma

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.

Input Data Blocks Needed by UFF


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 non-bond distance, the well depth of the Lennard-Jones potential, the scaling factor ζ, the effective charge, torsional barriers invoking a pair of sp3 atoms, torsional barriers involving a pair of sp2 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.

UFF Output Data Blocks


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=uffhessian0-0).

The file ufftopology

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 connectivity-block follows starting with the number of bond terms. Each line contains one bond term:

I   J     d   BO.

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:

J    I    K     wtyp    θ    nrJI   nrIK.

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:

wtyp = 1

linear case

wtyp = 2

trigonal planar case

wtyp = 3

quadratic planar case

wtyp = 6

octahedral case

wtyp = 9

all other cases.

θ is the angle value in degree. nrJI and nrIK 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:

I    J    K    L    nrJK     ttyp    ϕ    θIJK     θJKL.

Here are I,J,K and L the atom numbers. nrJK is the number of the bond between J and K. “ttyp” is the torsion type:

ttyp = 1

J (sp3)–K (sp3)

ttyp = 11

like ttyp=1, but one or both atoms are in Group 16

ttyp = 2

J (sp2)–K (sp3) or vice versa

ttyp = 21

like ttyp=2, but one or both atoms are in Group 16

ttyp = 22

like ttyp=2, but J or K is next a sp2 atom

ttyp = 3

J (sp2)–K (sp2)

ttyp = 9

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 NH3). In each line is one inversion term:

I    J    K     L    ityp1    ityp2    ityp3    ω1    ω2    ω3.

I,J,K and L are the atom numbers. Atom I is the central one. ityp1, ityp2, ityp3 are the types of the inversions:

ityp = 10

atom I is C and atom L is O

ityp = 11

like ityp=10, but L is any atom

ityp = 2

I is P

ityp = 3

I is As

ityp = 4

I is Sb

ityp = 5

I is Bi

ityp = 9

all other cases.

ω12 and ω3 are the values of the inversion angles in degree.

The nonbond terms follow starting with the number of the non-bonded terms. In each line is one nonbond term:

I   J     d .

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

        (      - γx)
˜λi = λi ⋅α + β ⋅e

and calculates a new Hessian with these modified eigenvalues.

20.2.6 Keywords for woelfling

Module WOELFLING reads options from data group


The below values of the options are default values with the following meaning:

 ninter           14

Number of interpolated structures for optimization.

 ncoord            2

Number of input structures provided by user.

 align             0

Align input structures by translation/rotation 0=yes, 1=no.

 maxit          40

Maximum number of iterations.

 dlst   3.00000000000000

Threshold for accuracy of LST-interpolation.

 thr  1.000000000000000E-004

Threshold for mean of norms of projected gradients.

 method q

Use standard optimization from initial LST-path (method q) or grow reaction path (method qg).


 riter 0

counts the number of completed iteration (no option).

20.2.7 Keywords for Modules dscf and ridft

$denconv real

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.d-7 or less.

$dft options

Listing of all possible sub-keywords for $dft (cross-references are given).

The user normally has to choose only the functional and the grid size, see below. All other parameters have proven defaults.

functional name

Specification of the functional, default is BP86, printed as functional b-p. 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):

   functional b-p

gridsize integer or minteger

Specification of the spherical grid (see section 20.2.7 on page 765). Default is gridsize m3.


   gridsize m3

gridtype integer

—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):

   gridtype 3

debug integer

Flag for debugging. debug 0 means no debug output (default), debug 1 means some output, debug 2 means a lot more output. Be careful!

nkk integer

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.


   nkk 3

ntheta integer
nphi integer

—not recommended for use—
Only for user-specified 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 user-defined Lobatto grids is the number of grid points, which must not exceed 2000 grid points.


  gridsize 9  
  ntheta 30  
  nphi 60


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.



radsize integer

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,

number of radial grid points = ioffrad+ (radsize- 1)* 5.

ioffrad is atom-dependent, 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).

diffuse integer

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

  gridsize m4  
  diffuse 2

rhostart integer
rhostop integer

—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,

   rhostart 50  
   rhostop 200

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:

gridsize 7  
radsize 14  
fullshell 1  
dgrenze 16  
fgrenze 16  
qgrenze 16  
fcut 16

Please refer to the corresponding sub-keywords for explanation.

If you want to use the reference grid, you have to skip the keyword gridsize, and type instead reference. Example:

   functional b-p  


Check if the selected grid is accurate enough for the employed basis-set by performing a numerical integration of the norm of all (occupied and virtual) orbitals. Useful for LHF. 771.

batchsize integer

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:

   functional b-p  
   gridsize m4  

symblock1 real
symblock2 real

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

xparameter integer

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

weight derivatives

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:

   gridordering 0

s-junc integer
p-junc integer

Within the semi-numerical integration scheme employed for local hybrid functionals, S- and P-junctions feature pre-screenings of the analytical two-center 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. s-junc and p-junc 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 P-junctions also depend on the grid weights. If larger grids are employed seeking more accurate results, the threshold for P-junctions should thus be adjusted (increasing the value to p-junc 6). The defaults settings are s-junc 6 and p-junc 6 for grad, rdgrad, and ridft. Too small values for s-junc and p-junc might result in non-converging calculations or deteriorated results. This holds particularly for excited state calculations. For escf the default values are therefore s-junc 8 and p-junc 8.

$electrostatic field

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.


$electrostatic field  
       0.1000E-03   0.000       0.000

$fermi tmstrt=<300.0> tmend=<100.0> tmfac=<0.9> hlcrt=<1.0E-01> stop=<1.0E-03> nue=<N>

Requests calculation of occupation numbers at a finite temperature T. For an orbital with the energy εi the occupation number ni [0,1] is calculated as

         (      )
n =  1erfc  εi --μ ,
 i   2      fT

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 HOMO-LUMO 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: 110-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 first-order SCF-calculation, i.e. perform only one SCF-iteration with the start MOs (which should be the orthogonalized MOs of two independent subsystems as is explained in detail in Chapter 18).

$fldopt options

Specification of options related with external electrostatic fields. The following options are available:

1st derivative on/off

Calculate numerical 1st derivative of SCF energy with respect to electrostatic field (default: off), increment for numerical differentiation is edelt (see below).

2nd derivative on/off

Calculate numerical 2nd derivative of SCF energy with respect to electrostatic field (default: off), increment for numerical differentiation is edelt.

edelt= real

Increment for numerical differentiation (default: 0.005).

fields on/off

Calculate SCF energy for non-zero external electrostatic
fields defined in $electrostatic field.

geofield on/off

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 non-negligible errors!!

$incore integer

By using this option the two-electron 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 two-electron integrals, i.e. SCF and DFT (including hybrid functionals); RHF and UHF.

The following condition must be met:

$scfdenapproxl 1

and rhfshells 1 or 2. It is advisable to set $thize as small as possible (e.g. $thize 0.1d-08) and to remove the keyword $scfdump.

Note: this keyword does not work for parallel runs.

$mo-diagram only nirreps=integer

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 mo-diagrams. 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 self-consistently. [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.

$mo output format format

If this line is present, the dscf program is forced to output the MOs using the new FORTRAN format format regardless of the format-option in data group $scfmo. Otherwise the input format will be used.

Example: $mo output format(3(2x,d15.8))

$natural orbitals

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

$natural orbital occupation

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. 

 a       1-5                     (     2.00000000000000 )  
 a       6                       (     1.99949836999366 )  
 a       7                       (     1.99687490286069 )  
 a       8                       (     1.00000000000000 )  
 a       9                       (      .00312509713931 )  
 a       10                      (      .00050163000634 )


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.

$restart dscf twoint

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 two-electron integrals from the files specified in $scfintunit, so there will be almost no loss of cpu-time.

All this information is normally provided by the previous dscf run if the keyword $scfdump (see there) was given.

$restartd data

Data provided by a previous dscf run that has been interrupted. This keyword is created when $scfdump was given.

$rundimensions data

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.


$scfconv integer

SCF convergency criterion will be 10-integer for the energy. Gradients will only be evaluated if integer > 6.

$scfdamp start=<.500> step=<.050> min=<.100>

Damping parameters for SCF iterations in order to reduce oscillations. The old Fock-operator 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 closed-shell 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.

$scfdebug options

Flags for debugging purposes. Following options are available:

vectors integer

Output level concerning molecular orbitals. integer=0 (default) means minimal output, >1 will output all start MOs and all MOs in each iteration.

density integer

Output level concerning difference density matrices.

debug integer

integer > 0 will dump a lot of information—be careful!

$scfdenapproxl integer

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.

$scfdiis options

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)

char=’FDS’ or ’SDF’ or ’FDS-SDF’

uses (FDS - SDF) as error vector.

char= none


char= sFDs

use S-12FDS12 - 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=FDS-SDF, maxiter=5, qscal=1.2, thrd=0.0, this implies DIIS damping in all iterations, prediag is switched of.


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

$scfintunit options

Disc space specification for two-electron 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:

 unit=30       size=35       file=twoint1  
 unit=31       size=35       file=/users/work/twoint2

Maximal 30 files may be specified in this way.

$scfiterlimit integer

Maximum number of SCF iterations (default: 30).

$scfmo none file=char

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.

format(format string)

This specifies the FORTRAN format specification which was used for MO output. The standard format is (4d20.14). (See data group $mo output format.)

Your data group $scfmo could look like this after a successful TURBOMOLE run :

$scfmo  scfconv=7  format(3(1x,d19.13))  
1  a1   eigenvalue=-.524127   nsao=6  
.1234567890123d+01 -.1234567890123d+00  .1234567890123d-01  
.1234567890123d+01 -.1234567890123d+00  
3  a2   eigenvalue=-.234810  

$scforbitalorder on/off

Order SCF MOs with respect to their energies (default: on)

$scforbitalshift options

To assist convergence, either the energies of unoccupied MOs can be shifted to higher energies or, in open-shell cases, the energies of closed-shell MOs to lower energies. In general a large shift may help to get better convergence.

Options are:


Automatic virtual shell shift switched off.

automatic real

Automatic virtual shell shift switched on; the energies of virtual orbitals will be shifted if the HOMO-LUMO 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 open-shell 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:

a1  1,2,4-6  (-.34)  
b1  8        (+.3)

$scftol real

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 = 10-(scfconv+1)-
   3⋅#bf (#bf = number of basisfunctions).

$scratch files

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:

$scratch files  
    dscf    dens          path1/file1  
    dscf    fock          path2/file2  
    dscf    dfock         path3/file3  
    dscf    ddens         path4/file4  
    dscf    statistics    path7/file7  
    dscf    errvec        path8/file8  
    dscf    oldfock       path9/file9  
    dscf    oneint        path10/file10

The first column specifies the program type (dscf stands for SCF energy calculations, i. e. the dscf program), the second column the scratch file needed by this program and the third column the pathname of the file to be used as scratch file.
$statistics options

The following options are allowed


Do not perform integrals statistics


Perform integrals statistics for dscf


see mpgrad

dscf parallel 


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 n4-step and may therefore be time-consuming) 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.

$thime integer

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)

$thize real

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.100E-04.


$closed shells

Specification of MO occupation for RHF, e.g. 

 a1g     1-4                                    ( 2 )  
 a2g     1                                      ( 2 )

$open shells type=1

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:

b2g     1                                      ( 1 )  
b3g     1                                      ( 1 )


This data group is necessary for ROHF calculations with more than one open shell. Example:

$rohf         1  
    a -a    a=0  b=0  
    h -h    a=1  b=2  
    a -h    a=1  b=2

This example is for the 7S state of chromium (3d54s1) in symmetry group I. Note that for this option being activated, $roothaan also has to be specified in your control file, although its parameter has no meaning in this case. For more details see Section 6.3.

For ROHF-calculations with only one open shell the Roothaan parameters a and b have to be specified within this data group (see also $rohf). Example:

    a = 3/4      b = 3/2

This example is for the 3P ground state of carbon (2p2) in symmetry group I. define recognizes most cases and suggests good Roothaan parameters.

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 open-shell electrons are simply ignored.


$alpha shells and $beta shells

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)

$alpha shells  
a      1-8                                    ( 1 )  
b      1-2                                    ( 1 )  
$beta shells  
a      1-7                                    ( 1 )  
b      1-3                                    ( 1 )


directs the program to carry out a UHF run. $uhf overwrites closed-shell occupation specification.

$uhfmo_alpha and $uhfmo_beta

These two data groups contain the UHF MO vectors for alpha and beta spin respectively (same syntax as $scfmo).


see $uhfmo_alpha

 functional b-p  
 gridsize m3

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 B-P86 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 so-called 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 b3-lyp_Gaussian and s-vwn_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 non-RI Hartree–Fock or DFT calculation.


Enforces an RI-J calculation if module ridft is used, can be used for Hartree-Fock as well as for DFT calculations with pure or hybrid functionals.


Obsolete keyword - use $rij instead!


Enforces a RI-JK calculation if module ridft is used, can be used for Hartree-Fock as well as for DFT calculations with pure or hybrid functionals.

$ricore integer

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)

$jbas file=auxbasis

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 RI-approximation for both Coulomb and HF-exchange (efficient for large basis sets). For this purpose needed (apart from $ricore):

$jkbas file=auxbasis

Cross reference for the file specifying the JK-auxiliary basis as referenced in $atoms. This group is created by the rijk menu in define.


Multipole-Accelerated-Resolution-of-Identity-J. This method partitions the Coulomb interactions in the near- and far-field parts. The calculation of the far-field part is performed by application of the multipole expansions and the near-field part is evaluated employing the RI-J 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.

  precision   1.0D-06  
  lmaxmom     10  
  nbinmax     8  
  wsindex     0.0  
  extmax      20.0  
  thrmom      1.0D-18

The following options are available:


specifies precision parameter for the multipole expansions. Low-precision MARI-J calculations require 110-6, which is the default. For higher precision calculations it should be set to 1 10-81 10-9.


maximum l-moment 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 far-field. 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.

Seminumeric HF-Exchange

If the keyword $senex is found in the control file, ridft performs a Hartree–Fock– SCF calculation using the seminumerical approximation for HF-exchange. Standard dft-grids can be used for the numerical integration. Smaller grids (-1,0) and the corresponding m-grids (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 dft-calculations. 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:

 gridsize 5


Use the Localized Hartree–Fock (LHF) method to obtain an effective Exact-Exchange Kohn–Sham potential (module dscf). The LHF method is a serial implementation for spin-restricted closed-shell and spin-unrestricted ground states.

 functional lhf  
 gridsize 6

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 test-integ to check if the selected grid is accurate enough for the employed basis-set.

How to do LHF runs

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)

   off-diag    on  
   numerical-slater off  
   pot-file save  
   asymptotic dynamic=1.d-3  
   homo     1b1u  
   homob    1b1u    # ONLY UNRESTRICTED  
   conj-grad conv=1.d-7 maxit=20 output=1 cgasy=1  
   slater-dtresh    1.d-9  
   slater-region     7.0 0.5 10.0 0.5  
   corrct-region             10.0 0.5  
   slater-b-region   7.0 0.5 10.0 0.5 # ONLY UNRESTRICTED  
   corrct-b-region           10.0 0.5 # ONLY UNRESTRICTED  
   correlation func=lyp

off-diag off

calculation of the KLI exchange potential. By default the LHF exchange potential is computed (off-diag on).

numerical-slater on

the Slater potential is calculated numerically everywhere: this is more accurate but much more expensive. When ECP are used, turn on this option.

numerical-slater off

leads to accurate results only for first-row elements or if an uncontracted basis set or a basis set with special additional contractions is used: in other cases numerical-slater on has to be used (this is default).


for asymptotic treatment there are three options:

asymptotic off

No asymptotic-treatment 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 density-matrix convergence or if Rydberg virtual orbitals are of no interest.

asymptotic on

Full asymptotic-treatment and use of the numerical Slater in the near asymptotic-region.

asymptotic dynamic=1.d-3

Automatic switching on (off) to the special asymptotic treatment if the differential density-matrix rms is below (above) 1.d-3. This is the default.

pot-file save

the converged Slater and correction potentials for all grid points are saved in the files slater.pot and corrct.pot, respectively. Using pot-file 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.

correlation func=functional

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 conjugate-gradient algorithm for the computation of the correction potential: rms-convergence (conj-grad conv=1.d-7), maximum number of iteration (maxit=20), output level output=0-3, asymptotic continuation in each iteration (cgasy=1).

With slater-dtresh= 1.d-9 (default) the calculations of the numerical integrals for the Slater potential is performed only if it changes more than 1.d-9.

Asymptotic regions specification:

corrct-region RFΔF
0RF - ΔF : basis-set correction potential
RF - ΔFRF + ΔF : smooth region
RF + ΔF + : asymptotic correction
Defaults: RF = 10, ΔF = 0.5

slater-region RNΔNRFΔF
0RN - ΔN : basis-set Slater potential
RN - ΔNRN + ΔN : smoothing region
RN + ΔNRF - ΔF : numerical Slater
RF - ΔFRF + ΔF : smoothing region
RF + ΔF + : asymptotic Slater
Note: RF - ΔF RF - ΔF
Defaults: RN = 7, ΔN = 0.5, RF = 10, ΔF = 0.5
Use correct-b-region and slater-b-region for the beta spin.

Two-component SCF (GHF)

Self-consistent two-component calculations (e.g. for spin-orbit interactions) can be carried out using the module ridft . The following keywords are valid:


enforces two-component-SCF calculations; this option is combinable with $rij, $rik and $dft.


switches on Kramers-restricted formalism


switches on collinear two-component formalism (not rotational invariant)


enforces DIIS for complex Fock operator.

All-electron relativistic approaches (X2C, BSS, DKH)

Relativistic all-electron 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.

$rdkh integer

switches on DKH calculation of order integer.

$dkhparam integer

selects parameterization of the DKH Hamiltonian. Valid values are 1 (=default), 2, 3, 4, and 5.

$dkhparam 1:
Optimum parametrization (OPT)
$dkhparam 2:
Exponential parametrization (EXP)
$dkhparam 3:
Square-root parametrization (SQR)
$dkhparam 4:
McWeeny parametrization (MCW)
$dkhparam 5:
Cayley parametrization (CAY)

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. Higher-order 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.

$rlocpara integer

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.

20.2.8 Keywords for Point Charge Embedding


Specification of position and magnitude of point charges to be included in the Hamiltonian. Each point charge is defined in the format

<x> <y> <z> <q>

with <x>, <y>, <z> being the coordinates and <q> its charge,e.g. 

$point_charges thr=<real> selfenergy nocheck list  
   2. 2. 2. 5.  
   5. 0. 0. 2.5

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 C1, 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)

20.2.9 Keywords for Periodic Electrostatic Embedded Cluster Method

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 three-dimensional system, 2 for a two-dimensional surface slab, and 1 for a one-dimensional 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.

  label x y z  

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.

  label x y z  

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.

  label charge  

Values of point charges (for each atom-type) , where label is the point charge label and charge specifies charge in atomic units.

  label charge  

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.


Well-separateness criterion for PFMM. Default is 3.0.


Minimum accuracy for lattice sums in PFMM. Default is 1.0d-8.

20.2.10 Keywords for COSMO

The Conductor-like 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:

  nppa= 1082  
  nspa=   92  
  disex= 10.0000  
  rsolv= 1.30  
  routf= 0.85  
  cavity closed  
  ampran= 0.1D-04  
  phsran= 0.0  
  refind= 1.3  
# the following options are not used by default  
  allocate_nps= 140  


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(ε) = ε-1
ε+x 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:

nppa= integer

number of basis grid points per atom
(allowed values: i = 10 × 3k × 4l + 2 = 12,32,42,92...)

nspa= integer

number of segments per atom
(allowed values: i = 10 × 3k × 4l + 2 = 12,32,42,92...)

disex= real

distance threshold for A matrix elements (┼ngstrom)

rsolv= real

distance to outer solvent sphere for cavity construction (┼ngstrom)

routf= real

factor for outer cavity construction in the outlying charge correction

cavity closed

  pave intersection seams with segments

cavity open

  leave untidy seams between atoms

ampran= real

amplitude of the cavity de-symmetrization

phsran= real

phase of the cavity de-symmetrization

refind= real

refractive index used for the calculation of vertical excitations and num. frequencies (the default 1.3 will be used if not set explicitly)


uses A matrix setup of TURBOMOLE 5.7


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:

# radii in Angstrom units  
o  1                                                   \  
   radius=  1.7200  
h  2-3                                                 \  
   radius=  1.3000

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

    cosmo      :  -0.003925  
    correction :   0.003644  
    total      :  -0.000282  
  ENERGIES [a.u.]:  
    Total energy            =      -76.0296831863  
    Total energy + OC corr. =      -76.0297567835  
    Dielectric energy       =       -0.0118029468  
    Diel. energy + OC corr. =       -0.0118765440  
    The following value is included for downward compatibility  
    Total energy corrected  =      -76.0297199849

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

DCCOSMO-RS: The DCOSMO-RS 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.

$dcosmo_rs file=filename.pot

activates the DCOSMO-RS method. The file defined in this option contains the DCOSMO-RS σ-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 25C are implemented in the current TURBOMOLE distribution (see parameter subdirectory):

































The DCOSMO-RS energies and total charges are listed in the COSMO section of the output:

    cosmo      :  -0.012321  
    correction :   0.011808  
    total      :  -0.000513  
 (correction on the COSMO level)  
  ENERGIES [a.u.]:  
    Total energy            =      -76.4841708454  
    Outlying charge corr. (COSMO)    =    -0.0006542315  
    Outlying charge corr. (DCOSMO-RS)=    -0.0011042856  
    Combinatorial contribution of the solute =    -0.0017627889  
    (at inf. dil. in the mixture/pure solvent. Not included in the total energy above)

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 DCOSMO-RS dielectric energy calculated form the corrected and the uncorrected COSMO charges, respectively (DCOSMO-RS). The charges are corrected on the COSMO level only. The Total energy includes the Ediel,RS defined in section 17.2. Additionally the combinatorial contribution at infinute dilution of the COSMO-RS 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 COSMO-RS contribution of the DCOSMO-RS energy depends on the reference state and the COSMO-RS parameterization (used in the calculation of the chosen COSMO-RS potential). Therefore, the DCOSMO-RS energies should not be used in a comparision with the gas phase energy, i.e. the calculation of solvation energies.

20.2.11 Keywords for Module riper

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:

$periodic n

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 (three-dimensional) lattice vectors need to be provided.

For 2D periodic systems two (two-dimensional) 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 n1 n2 n3

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 k1, k2 and k3 along the reciprocal lattice vectors b1, b2 and b3 as

k = k1b1 + k2b2 + k3b3.

For 2D periodic systems k3 = 0. In case of 1D periodicity k3 = 0 and k2 = 0. The three components kj (j = 1,2,3) of k are given as

kj = -i with i = - nj --1,- nj --1-+ 1,..., nj-- 1-- 1, nj --1.
     nj            2       2            2         2

with nj (j = 1,2,3) as integer numbers. For 3D periodic systems n1, n2 and n3 have to be specified. The component n3 can be omitted for 2D periodic systems. In case of 1D periodicity only the n1 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.

 kptlines 4  
 recipr 0.500 0.500 0.500 0.000 0.000 0.000 40  
 recipr 0.000 0.000 0.000 0.500 0.500 0.000 40  
 recipr 0.500 0.500 0.000 0.746 0.373 0.373 40  
 recipr 0.746 0.373 0.373 0.000 0.000 0.000 40

The calculated band structure is written to the file bands.xyz. Each line of the file contains five real numbers: the coordinates k1, k2 and k3 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:

  # general keywords  
   thrints    1.0d-12  
   lenonly    off  
   lchgprj    on  (for periodic systems)  
   lchgprj    off (for molecular systems)  
   northol    5  
   pqmatdiag  off  
   pqsingtol  1.0d-8  
  # CFMM control options  
   lmaxmom    20  
   nctrgt     10 (for periodic systems)  
   nctrgt     1  (for molecular systems)  
   wsicl      3.0  
   epsbext    1.0d-9  
   locmult    on  (for periodic systems)  
   locmult    off (for molecular systems)  
   locmomadd  2  
  # LMIDF control options  
   lpcg       on  
   lcfmmpcg   on  
   lmxmompcg  20  
   pcgtol     1.0d-9  
   pcgtyp     sp  
  # Gaussian smearing options  
   sigma      0.0d0  
   desnue     0.0d0  

The following options are available:

# general keywords


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.

# CFMM control options


Maximum l-moment 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 well-separateness criterion [107]. Octree boxes with centers separated more than sum of their lengths times 0.5×wsicl are considered as well-separated. 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.

# LMIDF control options


If set to on the low-memory 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 matrix-vector 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 l-moment 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 lower-precision 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.

pcgtyp char

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.

# Gaussian smearing options


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:

$pointvalper fmt=plt  
     orbs 2  
     k 3 2 1 a 1 i  
     k 0 0 0 b 2 r  
     nimg 2 2 1  
     npts 100 100 100  
     eps 5.0  
     ngrdpbx 50  

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 non-orthogonal 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 non-orthogonal 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:

$dosper width=real emin=real emax=real scal=real npt=integer  

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

20.2.12 Keywords for Modules grad and rdgrad

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

in $dft.

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.

20.2.13 Keywords for Module aoforce

This module calculates analytically harmonic vibrational frequencies within the HF- or (RI)DFT-methods for closed-shell and spin-unrestricted open-shell-systems. 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 RIDFT-energy 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.

General keywords


is the keyword for non-default options of gradient and second derivative calculations. Possibilities in case of the module aoforce are:

frequency analysis only

analysis only

to read a complete Hessian from the input file $hessian and perform only the frequency analysis

analysis [only] intcoord [print printlevel]

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.

$maxcor 50

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.

$forceconv 7

sets the convergence criterion for the CPHF-equations to a residual norm of 1.0d-7. Normally the default value of 1.0d-5 already provides an accuracy of vibrational frequencies of 0.01 cm-1 with respect to the values obtained for the convergence limit.

$forceiterlimit 10

fixes the maximum number of Davidson iterations for the solution of the CPHF-equations 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 C1 symmetry not to use 3N - 6(5) symmetry adapted but all 3N cartesian nuclear displacement vectors. This option may lead to a moderate speed-up 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 rotation-contributions, 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. 

 3 2.001  
 5 13

The interpolation then takes place between the mass(es) specified in $atoms (or the default mass(es), if none specified) and the mass(es) in $isosub. Take care of symmetry equivalent atoms, otherwise symmetry analysis will fail. This feature can not be used in a lowest eigenvalue search (keyword $les).
$isopts 6

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:


CPHF-iteration is done only for distortions, that are IR active.


CPHF-iteration 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 guess-Hessian and macro-iterations to find the solution vector(s) for the lowest eigenvalue(s). In each macro-iteration the CPHF-equations 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. 

 a1 3
 a2 all
 b2 1

The convergence criterion of the Davidson iterations for the solution of the CPHF-equations as well as the maximal residual norm for the lowest Hessian eigenvalue in the macro-iteration are specified by $forceconv as explained above.

The maximum number of macro-iterations is specified by  $lesiterlimit x

with the default x=25. The maximum number of iterations for each solution of the CPHF-equations is again determined by $forceiterlimit as shown above.

The convergence of the macro-iterations is strongly influenced by the size of the starting search-subspace. Generally all guess-Hessian eigenvectors corresponding to imaginary frequencies and at least two real ones are used as starting search-subspace. However it proved to be necessary to use even more vectors in the case of guess-Hessians with very large conditioning numbers.
$hesscond 8.0d-5
means that all eigenvalues with the quotient (eigenvalue)(max. eigenvalue) lower than 0.00008 are added to the starting search-subspace. Default is 1.0d-4.


Triggers the generation of input files for hotFCHT (program to calculate Franck-Condon-factors by R. Berger and co-workers). 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.

20.2.14 Keywords for Module evib

$dfdxi textout

can be used to generate text output of the matrix elements of the derivative of the Fock-operator. For bigger systems this can however generate very large output files. See 15

20.2.15 Keywords for Module escf

TDHF and TDDFT calculations

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 well-converged orbitals, specifying $scfconv 7 and $denconv 1d-7 for the ground-state 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 time-dependent DFT singlet or triplet or spin-unrestricted excitation energies (HF+RI(DFT))

$scfinstab rpas


$scfinstab rpat


$scfinstab urpa

2. TDA (for HF: CI singles) singlet or triplet or spin-unrestricted or spin-flip excitation energies (HF+RI(DFT))

$scfinstab ciss


$scfinstab cist


$scfinstab ucis


$scfinstab spinflip

3. Two-component TDDFT/TDA excitation energies of Kramers-restricted closed-shell systems

$scfinstab soghf


$scfinstab tdasoghf

4. Eigenvalues of singlet or triplet or non-real stability matrices (HF+RI(DFT), RHF)

$scfinstab singlet


$scfinstab triplet


$scfinstab non-real

5. Static polarizability and rotatory dispersion tensors (HF+(RI)DFT, RHF+UHF)

$scfinstab polly

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

$scfinstab dynpol nm  
  400 i

The number and symmetry labels of the excited states to be calculated is controlled by the data group $soes. Example:

b1g  17  
eu   23  
t2g  all

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 system-independent formatted logfile will be constructed for each IRREP. It can be re-used 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.

general keywords  

$rpacor n

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.

$spectrum unit

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. 

$cdspectrum unit

The calculated excitation energies and corresponding rotatory strengths are appended to a file named ’cdspectrum’. unit can have the same values as in $spectrum.

$start vector generation e

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

$velocity gauge

Enables calculation of dipole polarizability/rotatory dispersion in the velocity gauge. Active only for pure DFT (no HF exchange).

$sum rules unit
list of frequencies

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.

$rpaconv n

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.

$escfiterlimit n

Sets the upper limit for the number of Davidson Iterations to n. Default is n = 25.

GW Keywords

The main keyword that switches on a GW calculation. Provided that the response function is calculated setting this keyword will perform a standard G0W0 calculation with default values for the other flags. There are several options which can be added to the $gw keyword, the syntax is:

  nl 1  
  gam 0.001  
  output qpe.dat  

With the optional entries:

  nl <integer>

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.

  gam <real>

Default: 0.001. Infinitesimal complex energy shift. Negative value switches to calculating at that value but extrapolating to 0 in linear approximation.

  output <filename>

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.

BSE Keywords

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 quasi-particles energies.

  file <filename>

Default: qpenergies.dat. This option defines the file name of GW quasi-particle energies which are read in.


Default: Not set. Compute the auxiliary matrix for the static screened interaction iteratively instead of by Cholesky decomposition.

  thrconv <real>

Default: 1.0d-12. Threshold for convergence in case of the iterative computation of the static screened interaction.

  iterlim <integer>

Default: 100. Maximum number of iterations in case of the iterative computation of the static screened interaction.

20.2.16 Keywords for Module rirpa

The keyword $rirpa allows to specify the following options,

npoints n

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 RI-RPA.


Computes gradients in the direct RI-MP2 limit.

niapblocks n

Manual setting of the number of integral blocks, n, in subroutine rirhs.f. This is for developers; the default is determined with the $maxcor.

20.2.17 Keywords for Module egrad

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

$exopt n

which explicitly enforces that n-th excited state is optimized. n must not be larger than the number of states specified in $soes.


flag to compute Cartesian non-adiabatic 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 C1 symmetry.

20.2.18 Keywords for Module mpgrad

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

2) energy convergence $scfconv 6

$maxcor n

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.

 a1g     1-2  
 t1u     1

The data group $freeze specifies frozen orbitals, in the above syntax by irreducible representations. The symmetry-independent and for standard-applications recommended syntax is

 implicit core=5 virt=2

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:

$traloop n

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.

 type=intermed unit=61       size=0        file=halfint  
 type=1111     unit=62       size=0        file=moint#0  
 type=1112     unit=63       size=0        file=moint#1  
 type=1122     unit=64       size=0        file=moint#j  
 type=1212     unit=65       size=0        file=moint#k  
 type=1212a    unit=70       size=0        file=moint#a  
 type=gamma#1  unit=71       size=0        file=gamma#1  
 type=gamma#2  unit=72       size=0        file=gamma#2  
 type=1212u    unit=73       size=0        file=moint#u  
 type=1112u    unit=74       size=0        file=moint#v  
 type=gamma#1u unit=75       size=0        file=gamma#1u

The data group $mointunit specifies:

$statistics mpgrad

statistics run (estimation of disc space needed) for the adjustment of the file sizes will be performed.


calculation of MP2 pair correlation energies.

20.2.19 Keywords for Module ricc2

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 4-index integrals (for relaxed properties and gradients) on the settings for integrals storage in semi-direct SCF runs (i.e. $thime, $thize, $scfintunit). For the explanation of these keywords see Section 20.2.7.

$cbas file=auxbasis

Auxiliary basis set for RI approximation.


Freeze orbitals in the calculation of correlation and excitation energies. For details see Section 20.2.18.

$printlevel 1

Print level. The default value is 1.

$tmpdir /work/thisjob

Specify a directory for large intermediate files (typically three-index coulomb integrals and similar intermediates), which is different from the directory where the ricc2 program is started.

$maxcor n unit reference

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.

$spectrum unit

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. 

$cdspectrum unit

The calculated excitation energies and corresponding rotatory strengths are appended to a file named ’cdspectrum’. unit can have the same values as in $spectrum.


conv = 5

The purpose of this data group is twofold: It activates the Laplace-transformed implementation of SOS-MP2 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 Laplace-transformation.


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


mp2     d1diag  
cis(d)  energy only  
conv    = 8  
oconv   = 7  
lindep  = 15  
maxiter = 25  
mxdiis  = 10  
maxred  = 100  
iprint  = 1  
fmtprop = f15.8  
geoopt model=cc2 state=(a" 2)  
scs  cos=1.2d0   css=0.3333d0  

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

mp2 d1diag

Request the calculation of the D1 diagnostic in MP2 energy calculations (ignored in MP2 gradient calculations). Note that the evaluation of the D1 diagnostic increases the computational costs of the RI-MP2 energy calculation roughly by a factor of 3.

cis(d) energy only

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 CCR0--1--1---0 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 coupled-cluster 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 =conv-1, 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 one-electron 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 open-shell 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 D1 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 D1 diagnostic is evaluated by default (without significant extra costs).


calculates the second-order corrections to the CCSD(T) energy from the interference-corrected MP2-F12 (INT-MP2-F12) 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)-INT-F12 energy is printed. The intcorr all keyword writes on the output all pair energies.



ansatz char

char=1, 2* or 2
The ansatz flag determines which ansatz is used to calculate the RI-MP2-F12 ground state energy.
(Ansatz 2 is used if ansatz is absent.)

r12model char

char=A, A’ or B
The r12model flag determines which approximation model is used to calculate the RI-MP2-F12 ground state energy.
(Ansatz B is used if r12model is absent.)

comaprox char

char=F+K or T+V
The comaprox flag determines the method used to approximate the commutator integrals [T,f12].
(Approximation T+V is used if comaprox is absent.)

cabs char val

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.0d-08 is used if cabs is absent.)

examp char

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 orbital-invariant method. For fixed and noinv only the diagonal amplitudes are non-zero 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 open-shell calculations noflip supresses the use of spin-flipped geminal functions.
(The fixed flip method is used if examp is absent.)

pairenergy char

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 RI-MP2-F12 correlation energies is always printed out.

corrfac char

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.

cabsingles char

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

r12orb char

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 ROHF-based UMP2 calculations rohf orbitals can be used, which also implies that the $freeze data group options refer to ROHF rather than semi-canonical orbitals. For closed-shell species, localised orbitals can be used with either the Boys or Pipek-Mezey 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 f12, f12g12, f122 and f12r12. 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.


irrep=au  multiplicity=1 nexc=4  npre=6  nstart=8  
irrep=bg  multiplicity=3 nexc=2  npre=4  nstart=5  
spectrum states=all operators=diplen,dipvel  
tmexc    istates=all fstates=all operators=diplen,dipvel  
exprop   states=all operators=qudlen  
twophoton states=all operators=(diplen,diplen) freq=0.1d0  
momdrv   states=all operators(diplen,soc) freq=0.0d0  
xgrad    states=(ag{3} 1)  
conv    = 6  
thrdiis = 2  
preopt  = 3  

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 two-photon transition moments between ground and excited states. It recognizes the optional suboptions states for specifying the states for which the two-photon transition moments should be computed, operators for specifying a list of pairs of one-electron transition operators, and freq for a list of frequencies for one of the photons. If these suboptions are specified, two-photon 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 two-photon transition moments should be computed, operators for specifying a list of pairs of one-electron transition operators. For phosphorescence lifetimes operators should be set to (diplen,soc) and freq to 0.0d0.


requests the calculation of first-order 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 non-linear 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 %T1 and %T2 diagnostics which were identical with those used in the CC response code of the Dalton program.


fop unrelaxed_only operators=diplen  
sop operators=(diplen,diplen) freq=0.077d0  
conv = 6  
zconv = 6  
thrsemi = 3

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 first-order properties (expectation values). The operators flag can be followed by a list of operators (see below) for which the first-order 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 orbital-relaxed first-order properties, which require solution the CPHF-like Z-vector equations. Default is the calculation of unrelaxed and orbital-relaxed first-order 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 second-order 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 symmetry-allowed elements of the dipole-dipole 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 orbital-relaxed approach. Note that the orbital-relaxed approach can not only be used in the static limit (freq=0.0d0). For further restrictions for the computation of second-order 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 closed-shell 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 semi-canonical formulation for the calculation of (transition) one-electron densities. Switched on by default. The semi-canonical formulation is usually computationally more efficient than the non-canonical formulation. Exceptions are systems with many nearly degenerate pairs of occupied orbitals, which have to be treated in a non-canonical way anyway. (See also explanation for thrsemi below).


use non-canonical formulation for the calculation of (transition) one-electron densities. Default is to use the semi-canonical 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 non-canonical fashion will be set to 10-thrsemi. If set to tight the semi-canonical algorithm will become inefficient, if the threshold is to large the algorithm will become numerical unstable


threshold for preoptimizating the so-called Z vector (i.e. the lagrangian multipliers for orbital coefficients) with a preceding RI-CPHF calculation with the cbas auxiliary basis. The RI-CPHF 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 RI-JK Hartree-Fock calculation.


disable the preoptimization of the Z vector by a preceding RI-CPHF calculation with the cbas basis set. (Note that the preoptimization is automatically deactivated if the ricc2 calculation is based on a foregoing RI-JK Hartree-Fock calculation.)

Common options for keywords in the data groups $ricc2, $response, and $excitations:


input of operator labels for first-order 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: μ|riO|ν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 μ|riOrjO|ν
(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) Mij, the isotropic second moment α = 1
3trM and the anisotropy

    ┌ --------------------------------
    ││   ∑z                   ∑z
β = ∘ 1   (Mii - Mi+1,i+1)2 + 3  M 2i,i+1.
      2 i=x                   i=x

Furthermore the traceless quadrupole moment

Θij = 1⟨3rirj - r2δij⟩

(including nuclear contributions) is given.


angular momentum μ|LiO|ν
(all three components, individual components can be specified with the labels xangmom, yangmom, zangmom).


electronic force on nuclei μ|ZIrIi
 rI3|ν, where ZI is the charge of the nucleus I and rI 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 first-order 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,3-5; b1u{1} 1-3; b2u4). This will select the triplet ag states no. 1, 3, 4, 5 and the singlet b1u states no. 1, 2, 3 and the singlet (which is default if no {} is found) b2u state no. 4.

istates=all fstates=all

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 double-substitution-based diagnostics D2.


Write MP2/CC2 natural occupation numbers and natural orbitals to a file.

$cgrad 1000

Calculate the error functional δRI for the RI approximation of (ai|bj) integrals

     1 |⟨ab||ij⟩exact - ⟨ab||ij⟩RI|2
δRI = 4----ϵ---ϵ-+-ϵ-- ϵ-----
            a   i  b   j

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 RI-MP2 and RI-CC2 calculations (see Section 1.5).

20.2.20 Keywords for Module ccsdf12

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:

ccsdapprox label

defines the approximation to CCSD-F12 which will be used if the MP2-F12 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]
CCSD-F12b ccsd-f12b
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 CCSD-F12b calculations also the CCSD-F12a 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.

20.2.21 Keywords for Module pnoccsd

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.

$cbas file=auxbasis

Auxiliary basis set for RI approximation.


Freeze orbitals in the calculation of correlation and excitation energies. For details see Section 20.2.18.

$printlevel 1

Print level. The default value is 1.

$tmpdir /work/thisjob

Specify a directory for large intermediate files (typically three-index coulomb integrals and similar intermediates), which is different from the directory where the program is started.

$maxcor n unit reference

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.


conv = 1

Only needed for test purposes. It sets the accuracy for the numerical Laplace-transformation to 10-conv. For OSV-PNO-MP2 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.


localize boys  
prepno   davidson  
mxrdim   800  
tolpno=   1.00E-7  
tolosv=   1.00E-8  
tolri=    1.21E-3  
tolpair=  4.64E-6  
opnos on  
tolcapno= 1.00E-9   3.16E-8  
tolosc=   1.00E-10  3.15E-9  
tolopno=  1.00E-9  
toloso=   1.00E-10  
projector 3  
conv    = 7  
oconv   = 3  
lindep  = 12  
maxiter = 25  
mxdiis  = 10  
maxred  = 100  
scs  cos=1.2d0   css=0.3333d0  


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 Foster-Boys, pm for Pipek-Mezey, ibo for intrinsic bond orbitals (IBOs) and none for canonical (deprecated, only meant for testing) orbitals. Default are Foster-Boys orbitals, IBOs are recommended if seperation of σ- and π-type orbitals is important. (Pipek-Mezey orbitals become very delocal with diffuse basis sets.)


switch between to algorithms for the generation of OSV coefficients. Possible choices are full for an O(N4) scaling direct diagonalization (cmp. [237]) and davidson for an O(N3) 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 PNO-MP2 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 PNO-MP2-F12 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 pair-specific auxiliary basis sets for the local RI approximation. If not given tolri is set to 10712 ×√tolpno-, which is the recommended value.


specifies the energy threshold for selecting the significant pairs. If not given tolpair is set to (0.1 ×tolpno)23


enables (on) or disables (off) for F12 calculations the use of OPNOs for the occupied orbital spaces in the projectors for the three-electron 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 three-electron 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 three-electron 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 AQij for the three-electron integrals. Possible choices are:
projector 1 for 1 -O1O2 -V 1V 2 - O1V 2′′-V 1′′O2
projector 2 for 1 -O1O2 -V 1V 2 -O1V 2′′- V 1′′O2
projector 3 for 1 -O1O2 -V 1V 2 -O1V 2′′-V 1′′O2
With projector 3 AQij becomes independend of the system size and is therefore the default. But note that for the O(N4) 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 pre-PNOs 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.

20.2.22 Keywords for Module relax

$optimize options

define what kind of nonlinear parameters are to be optimized by relax and specify some control variables for parameter update.

Available options are:

internal on/off

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

redundant on/off

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

cartesian on/off

optimize molecular structures in the space of (symmetry-distinct) cartesian coordinates (default: off).

basis on/off suboptions

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

global on/off

optimize a global scaling factor for all basis set exponents (default: off).

  • basis and global have to be used exclusively!
  • if $optimize has been specified but $forceapprox is absent, the option $forceinit on is switched on by default.
  • specification of the option $interconversion on will override $optimize!
$coordinateupdate options

define some variables controlling the update of coordinates.

Available options are:

dqmax real

maximum allowed total change for update of coordinates. The maximum change of individual coordinate will be limited to dqmax2 and the collective change dq will be damped by dqmaxdqdqif dqdq> dqmaxq
(default: 0.3)

interpolate on/off

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)

statistics on/integer/off

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

$gdiishistory file=char

the presence of this keyword forces relax to provide informational output about the usage of DIIS for the update of the molecular geometry.

$interconversion options default=off

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

on/off options

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:

cartesian –> internal coordinate gradient hessian

cartesian <– internal coordinate

the direction of the transformation is indicated by the direction of the arrow

Note: specification of $interconversion on will override $optimize!
$forceupdate method options

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)

ms suboptions

Murtagh–Sargent update

dfp suboptions

Davidon–Fletcher–Powell update

bfgs suboptions

Broyden–Fletcher–Goldfarb–Shanno update

dfp-bfgs suboptions

combined (bfgs+dfp) update

schlegel suboptions

Schlegel update

ahlrichs suboptions

Ahlrichs update (macro option)

suboptions if method=ms, dfp, bfgs, schlegel, ahlrichs


number of structures used


maximum number of geometries (= rank of the update procedure, for ahlrichs only)


minimum number of geometries needed to start update

if method=ms, dfp, bfgs: maxgeo=2, mingeo=1 as default

additional suboptions if method=ahlrichs

modus= char fmode

for an explanation see suboptions pulay given below e.g. ahlrichs numgeo=7 mingeo=3 maxgeo=4 modus=<g|dg> dynamic


if the macro option ahlrichs has been chosen and n=numgeo, ncycl=‘number of geometries available’
  • if ncycl < n: geometry update by inter/extrapolation using the last two geometries
  • if ncycl n: diagonal update for the hessian by least mean squares fit; pulay update for the geometry (using specified modus, fmode (see pulay below))
  • if (ncycl max(5,n + 3) and max(g) < 0.01 and g < 0.001 ) or Hij0ij: diagonal update is replaced by multidimensional BFGS (rank n) update for the hessian

pulay suboptions

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:


number of geometries to be utilized


maximum number of geometries


minimum number of geometries needed to start update

modus=char fmode

char=<g|g> or <g|dq> or <dq|dq> defines the quantity to be minimized (dq = internal coordinate change).
fmode specifies the force constants to be used (only if char=<g|dq> or <dq|dq>!)
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 |gg⋅|*dq|dq| > -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=<dq|dq> static fail=0.2

options for $forceupdate


update only the diagonal force constants (update for off-diagonals will be suppressed) (only active if method=ms, dfp, bfgs)

offdamp real

this allows to damp off-diagonal force constants by 1/real (compare offreset, which discards off-diagonals 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 non-diagonal update has lead to too large non-diagonal elements of the hessian.


reset off-diagonal 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 off-diagonals!


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

$forceinit option

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

  • internal coordinates (supplied in $intdef ... fdiag=..)
  • a global scale factor ( $global ... fdiag=..)

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 B-matrix with respect to cartesian coordinates, which is needed to transform H(cart) H(int) exactly by specifying no dbdx.

$last SCF energy change = real

$last MP2 energy change = real

These keywords depend on the optimization task to be processed and are updated by the corresponding program (i. g. SCF energy).

$m-matrix options

This data block contains non-default specifications for the m-matrix diagonals. This is of use if some cartesian atomic coordinates shall be kept fixed during optimization.

Available options are:

integer real real real

atomic index followed by diagonal elements of the m-matrix for this atom

$scratch files

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:

       $scratch files  
        relax   ftmp          path/file  

The first column specifies the program, the second column the scratch file and the third column the pathname of the file to be used as scratch file.

Input Data Blocks Needed by RELAX

$intdef or $redundant

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.

#  coordinate    increment  
     1            0.0600  
     8           -0.0850

$forceapprox options

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 10-digit f10.x formats (e.g. x=4,6,..) are possible and will be used, after being manually specified within $forceapprox. See the example below:

$forceapprox  format=(8f10.4)  
    -.0108    0.3347  
    0.2101    0.0299    1.3347  
    0.0076    0.1088    0.0778    0.6515

$hessian (projected)

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

RELAX Output Data Groups


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 quasi-Newton type geometry optimizations; this matrix will be improved in subsequent optimization cycles if one of the variable-metric 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 DIIS-type geometry optimizations. It will be initialized by relax specifying: $forceupdate pulay modus=<dq|dq> 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.

$maximum norm of cartesian gradient = real

$maximum norm of internal gradient = real

$maximum norm of basis set gradient = real

real is the absolute value of the maximum component of the corresponding gradient.

Other Input/Output data used by RELAX

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.


20.2.23 Keywords for Module statpt

  itrvec      0  
  update      bfgs  
  hssfreq     0  
  hssidiag    0.5  
  radmax      0.3  
  radmin      1.0d-4  
  tradius     0.3  
  threchange  1.0d-6  
  thrmaxdispl 1.0d-3  
  thrmaxgrad  1.0d-3  
  thrrmsdispl 5.0d-4  
  thrrmsgrad  5.0d-4

Only non-default values are written in the control file except:

  itrvec 0

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.0d-4).


Initial value for trust radius (default tradius: radmax = 0.3).

Convergence criteria  


threshold for energy change (default: 1.0d-6).


threshold for maximal displacement element (default: 1.0d-3).


threshold for maximal gradient element (default: 1.0d-3).


threshold for RMS of displacement (RMS = root mean square)
  (default = 5.0d-4)


threshold for RMS of gradient (default: 5.0d-4).

All values are in atomic units.

20.2.24 Keywords for Module moloch

$properties specifies the global tasks for program moloch by virtue of the following options

    trace                              off  
    moments                            active  
    potential                          off  
    cowan-griffin                      off  
    localization                       off  
    population analyses                off  
    plot                               off  
    firstorder                         off  
    fit                                off

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 cowan-griffin), there are additional data groups allowing for more detailed specifications, as explained below.


if moment is active you need

  0th 1st 2nd 3rd  
  point .0 .0 .0

to compute the 0th, 1st, 2nd and 3rd moment at the reference point 0 0 0.


if potential is active you need

$points #1  
  pot fld fldgrd shld  
  point .0 .0 .0

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 boys-localization 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. Non-defaults may be specified using the following suboptions:

lmofile= filename

locxyz dir1 dir2 dir3

threshold real

sweeps integer

population analyses

if population analyses is active you need

 spdf molap netto irpspd irpmol mommul

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


to perform a L÷wdin population analysis (options are invalid here). A L÷wdin population analysis is based on decomposing √S--D√S- instead of DS in case of a Mulliken PA.


 momao maodump maofile=mao all

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

maofile=mao all

print all MAOs to file mao.

This kind of population analysis basically aims at so-called shared electron numbers (SEN) between two or more atoms. By default 2-, 3- and 4-center 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
 2-center  threshold =  real
 3-center  threshold =  real
 4-center  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,8-10,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).


$mao selection  threshold= 0.09  
   atom c 1,3-5  nmao= 5  method= eig  threshold= 0.1  
   atom o 2      nmao= 3  method= man olabel  
 olabel  3-5


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

$grid     #1  
 mo 4a1g  
  origin      .000000      .000000      .000000  
 vector1     1.000000      .000000      .000000  
 vector2      .000000     1.000000      .000000  
 grid1 range    -5.000000     5.000000 points   100  
 grid2 range    -5.000000     5.000000 points   100  
 outfile = 4a1g

to obtain two-dimensional 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 old-fashioned way to plot MOs and densities. A new—and easier—one is to use $pointval, as described below.


if fit is active you need

shell     number_of_gridpoints     distance_from_vdW_surface  
refine    value_of_potential


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.

element_symbol    van_d_waals_radius

One line per element has to be specified, it contains the name of the element and the van der Waals radius in [Bohr].

20.2.25 Keywords for wave function analysis and generation of plotting data

Properties of RHF, UHF and (two-component) GHF wave functions as well as those of SCF+MP2 densities or such from excited state DFT-calculations can be directly analyzed within the respective programs (dscf, ridft, mpgrad, rimp2 and egrad). In case of spin-unrestricted 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 MP2-corrected 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 < p2 >,< p4 > and the Darwin term ( 1∕ZA * ρ(RA)).


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:

$moments <i>  
x1 y1 z1  
x2 y2 z2  

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 spin-unrestricted 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.

Two-component wavefunctions (only module ridft and only if $soghf is set): In two-component calculations instead of Sz |(Sx,Sy,Sz)| is written to the output. Additionally a vector-file 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 px,py,pz (and analogous for d and f functions) is displayed (lengthy output).

atoms list of atoms

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.

mosum list of MOs

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.

mo list of MOs

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

dos width=real points=integer

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_a-b, dos_alpha, dos_beta; for UHF cases) contain energies (first column), resulting DOS for the respective energy (second column) as well as s-, p-, d-contributions for the respective energy (following columns).


$pop mo 23-33 dos atoms 2,3,7-8

leads to Mulliken PA (CAO-basis) for each of the eleven MOs 23-33, regarding only contributions from atoms 2-3 and 7-8 (results are written to standard output) and generation of file(s) with the respective simulated density of states.

$pop nbo

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 tw to circumvent numerical difficulties in computing Ow

(default: tw=1.d-6).

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.


$pop nbo AO ab short atoms 1,2,6

leads to a natural population analysis (AO-basis) 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 $nbonmb-block in the control file. Example:

  ni s:4 p:2 d:1  
  o  s:2 p:1

leads to a NMB set for Ni of 4 s-, 2 p- and 1d-functions and for O of 2 s- and 1 p-functions.

$pop paboon

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+MP2-density in case of MP2 and the GHF total density for (two-component-)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.

$mao selection options

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 B-Ne, 6 for Na-Mg, etc. Exceptions are Elements Sc (Y, La), Ti (Zr, Hf), V (Nb, Ta) for which not all five d-shells are included, but only 2, 3 or 4, respectively. This modification leads to better agreement with partial charges calculated by an ESP-fit.

thr <real>

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 ESP-fit. 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,3-4 eig 7

leads to choice of the 7 MAOs with largest eigenvalue at atoms 1, 3-4.


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

mo list of MOs

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

pipmez or pm or pipek-mezey

Pipek-Mezey localization is used.


Foster-Boys localization is used.


Intrinsic Bond Orbitals are used.


Compute LMO centers, i.e. the diagonal matrix elements of the position operator ⃗ri = i|⃗r|i


Compute LMO spreads defined as σi = ∘ ⟨i|(⃗r---⃗ri)2|i⟩ = ∘ ⟨i|⃗r2|i⟩2 --⃗r2
            i, where ⃗ri 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 single-point 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 two-component 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, Bragg-Slater radii, rBS, are taken as shell radii, for each atom the number of points is given by 1000 rBS2, the total number of points is the sum of points for each atom reduced by the number of points of overlapping spheres. Non-default shells (one or more) can be specified as follows:

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), 129-145 (1984)) may be obtained by

$esp_fit kollman

Here five shells are placed around each atom with r=1.4*rvdW + k, k=0pm, 20pm, 40pm, 60pm, 80pm, and rvdW are the van-der-Waals radii of the atoms.


drives the calculation of space-dependent molecular quantities at 3D grids, planes, lines or single points. Without further specifications the values of densities are plotted on a three-dimensional 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 RHF-dscf/ridft calculations you get the total density on file td.plt, for UHF-dscf/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

$pointval integrate x y z r eps

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

mo list of MO numbers

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

Two-component 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 spinor-parts are calculated, whose weights (the probability of finding the electron in this part) lie above the threshold.

lmo list of LMO numbers

calculation of amplitudes of LMOs (previously generated by $localize) ordered by the corresponding diagonal element of the Fock matrix in the LMO basis.

nmo list of NMO numbers

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 Kohn-Sham exchange-correlation potential. It is only valid for DFT calculations and it works for all exchange-correlation functionals, including LHF. Note that for hybrid functionals, only the Kohn-Sham part of the potential will be computed (the HF part is a non-local-operator and can’t be plotted). For GGA functional the full potential will be computed (local and non-local terms)

For line plots the output file is tx.vec . For UHF calculations the output files are tx.vec (alpha-spin potential) and sx.vec (beta-spin potential).

For a line plot the file has three columns: 1: total potential 2: local term ( or Slater-potential for LHF) 3: non-local 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 (E-field) 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 (E-field) 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,vx,vy,vz.


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. Self-explaining 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 α ⃗v (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 two-dimensional 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 non-default boundarys, basis vector directions, origin and resolutions may be specified as follows:

grid1 vector 0 3 0 range -2,2 points 200  
grid2 vector 0 0 -7 range -1,4 points 300  
grid3 vector 1 0 0 range -1,1 points 300  
origin 1 1 1

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 4-5 and 9 in the following example:

  grid1 vector 0 3 0 range -2,2 points 200  
  grid2 vector 0 0 -7 range -1,4 points 300  
  orient 4-5,9  
  rotate 25 3

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:

$pointval geo=plane  
grid1 vector 0 1 0 range -2,2 points 200  
grid2 vector 0 0 1 range -1,4 points 300  
origin 1 1 1

Values are calculated at a plane spanned by vectors (0,1,0) and (0,0,1) centered at (1,1,1).

$pointval geo=line  
grid1 vector 0 1 0 range -2,2 points 50  
origin 0 0 1

Values are calculated at a line in direction (0,1,0) centered at (0,0,1). Output format as above.

$pointval geo=point  
7 5 3  
0 0 7

Values are calculated at the two points (7.0,5.0,3.0) and (0.0,0.0,7.0).

Plane-averaged density can be computed by

$pointval dens averagez fmt=vec  
grid1 vector 1 0 0 range -10,10 points 100  
grid2 vector 0 1 0 range -10,10 points 100  
grid3 vector 0 0 1 range -20,20 points 200  
origin 0 0 0

The generated file td.vec will contain the quantity

      ∫ ∫
ρ(z) =     dxdyρ(x,y,z)

20.2.26 Keywords for Module frog

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 online-help 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:

$nsteps 75

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.

$natoms 9

Number of atoms in system.

$current file=mdlog.aa

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.

$log file=mdlog.ZZ

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

$turbomole file=control

Where to look for TURBOMOLE keywords $grad etc.


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:

  canonical T=500 t=100  
  from t= -25.0000000000       until t=  0.00000000000

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 2-10 times larger than the time step. Note that user-defined actions are presently not supported in canonical dynamics mode.

These are optional keywords:

$seed -123

Integer random number seed


Arbitrary title

 100                mdlog.P  
 71                 mdlog.Q

  length        50  
  response      1

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 angstroms  
  type        elps  
  limits      5.0 10.0 7.5  
  constant    2.0  
  springlen   1.0  
  temperature 300.0

$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, semi-major 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 quasi-temperature 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.

$constraints angstroms  
  tolerance   0.05  
  adjpercyc   0.25  
  type H O 0.9 1.2  
  type F C 0.0 1.7  
  type H C -1.0 1.2  
  2 1 0.0  
  3 1 1.54  
  4 1 -1.0


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.

type X A const rmax

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 H-C distance multiplied by the absolute value of const.

Explicitly specified constraints are listed by atom index and supercede auto-generated 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.

User-defined 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:

     fix total energy from t=2000.0  
     anneal from t=2500.0  
     free from t=3000.0

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:

     fix temperature from t=<real>  
     fix total energy from t=<real>

These commands cause velocities to be scaled so as to keep the average kinetic energy (i.e. quasi-temperature) 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).

     set temperature at t=<real> to x=<real> K  
     set total energy at t=<real> to x=<real> H  
     set kinetic energy at t=<real> to x=<real> H  
     set position file=<filename> at t=<real>  
     set velocity file=<filename> at t=<real>  
     set velocity at t=<real> random  
     set velocity at t=<real> zero

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.

     anneal from t=<real>  
     anneal from t=<real> x=<real>  
     quench from t=<real>  
     quench from t=<real> x=<real> file=<file>  
     relax at t=<real>

In Simulated Annealing MD, the temperature of a run is lowered so as to find minimum-energy 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 non-dynamical geometry optimization can be carried out in a subdirectory (relax).

     free from t=<real>

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 Tully-type 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

 fix total energy    from t=  0.00000000000

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 non-adiabatic couplings; this keyword requires the use of weight derivatives in section dft


keyword needed to collect Cartesian non-adiabatic coupling vectors along the trajectory

$exopt 1

keyword needed to ensure dynamics starting in S1

$ex_energies file=ex_energies

collects excitation energies along the trajectory

$integral_ex file=integral_ex

collects time integration of excitation energies along the trajectory

$sh_coeffs file=sh_coeffs

collects amplitudes of the adiabatic states along the trajectory

$nac_matrix file=nac_matrix

collects NAC elements 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 Tamm-Dancoff approximation (TDA) to TDDFT (see 8). In addition an optional keyword for the md_master file can be used:

$gap_threshold <real>

enforces a switch to the ground state in case the S1-S0 energy gap drops below <real> eV. As default a switch to S0 is enforced if the S1 TDDFT-TDA 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.

20.2.27 Keywords for Module mpshift

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.

$traloop n

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.

$csconv real

Sets the convergence threshold for the shielding constant of succeeding CPHF iterations. The unit is ppm and the default value is 0.01.

$csconvatom integer

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

$thime, $thize, $scftol, $scfintunit, $scfmo

have the save meaning as in dscf (see Section 20.2.7);
Since mpshift works ’semi-direct’ it uses the same integral storage.

$scratch files

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:

$scratch files  
    mpshift   csssmat         path1/file1  
    mpshift   cshsmat         path2/file2  
    mpshift   csdgsmat        path3/file3  
    mpshift   csusmat         path4/file4  
    mpshift   dens            path5/file5  
    mpshift   fock            path6/file6  
    mpshift   dfock           path7/file7  
    mpshift   idvds1          path8/file8  
    mpshift   idvds2          path9/file9  
    mpshift   idvds3          path10/file10  
    mpshift   jdvds1          path11/file11  
    mpshift   jdvds2          path12/file12  
    mpshift   jdvds3          path13/file13  
    mpshift   cshmmat         path14/file14

$trast , $trand traloop-number

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

    $trast <number>  
    $trand <number>

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.


UaiBβ will be written to file for a subsequent calculation of VCD rotational strengths by aoforce.

20.2.28 Keywords for Parallel Runs

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 V-Class, SP3-SMP and HP S/X-Class


for systems with fast communication like Fast-Ethernet, 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.

$scratch files  
   dscf  dens       /home/dfs/cd00/cd03_dens  
   dscf  fock       /home/dfs/cd00/cd03_fock  
   dscf  dfock      /home/dfs/cd00/cd03_dfock  
   dscf  ddens      /home/dfs/cd00/cd03_ddens  
   dscf  xsv        /home/dfs/cd00/cd03_xsv  
   dscf  pulay      /home/dfs/cd00/cd03_pulay  
   dscf  statistics /home/dfs/cd00/cd03_statistics  
   dscf  errvec     /home/dfs/cd00/cd03_errvec  
   dscf  oldfock    /home/dfs/cd00/cd03_oldfock  
   dscf  oneint     /home/dfs/cd00/cd03_oneint

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
$2e-ints_shell_statistics    file=DSCF-par-stat
$2e-ints’_shell_statistics    file=GRAD-par-stat

The statistics files have to be generated with a single node dscf or grad run. For a dscf statistics run one uses the keywords:

$statistics  dscf parallel  
$2e-ints_shell_statistics    file=DSCF-par-stat  

and for a grad statistics run:

$statistics  grad parallel  
$2e-ints’_shell_statistics    file=GRAD-par-stat  

maxtask is the maximum number of two-electron integral tasks,
maxdisk defines the maximum task size with respect to mass storage (MBytes) and
dynamic_fraction is the fraction of two-electron integral tasks which will be allocated dynamically.