21.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. They are sorted according to program. Complete control files are provided in Chapter 22. An alphabetical list of all keywords is given in the index.

21.2.1 General Keywords

$operating system unix  
$path  
$lock off  
$suspend off

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

$statistics dscf
or
$statistics mpgrad

Only a statistics run will be performed to determine file space requirements as specified for dscf or mpgrad. On return the statistics option will be changed to $statistics off.

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

$basis


specification of basis sets.

$ecp

specification of effective core potentials.

$jbas


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.

$forceapprox


approximate force constant for geometry optimizations.

The control file must end with this keyword:

$end

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

$title


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.

$atoms


Example:

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

RHF

$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


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 22.6 and the tables of Roothaan parameters in Section 6.3.

UHF

$uhf

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.

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

stre

bond stretch (default: 0.5)

invr

inverse bond stretch (default: 0.5)

bend

bond angle (default: 0.2)

outp

Out of plane angle (default: 0.2)

tors

dihedral or “torsional” angle (default: 0.2)

linc

Special angle coordinate for collinear chains, bending of the chain a–b–c in the plane of b–c–d (default: 0.2)

linp

bending of the chain a–b–c perpendicular to the plane of b–c–d
(default: 0.2)

wstr

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)

winv

inverse stretch of a weak bond (default: 0.05)

wbnd

bond angle involving at least one weak bond (default: 0.02)

wout

Out of plane angle for weak bonds (default: 0.02)

wtor

dihedral angle for weak bonds (default: 0.02)

wlnc

linc coordinate for weak bonds (default: 0.02)

wlnp

linp coordinate for weak bonds (default: 0.02)

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

maxcycle


number of max. optimization cycles (maxcycle=1: single-point calculation).

modus


can have the values +1 or -1. If modus = -1 only the topology will be calculated.

nqeq

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.

iterm


switch for the different types of force field terms:

100000

bond terms will be calculated.

010000

angle terms will be calculated.

001000

torsion terms will be calculated.

000100

inversion terms will be calculated.

000010

non bonded van der Waals terms will be calculated.

000001

non bonded electrostatic terms will be calculated.

econv, gconv


convergence criteria for energy and gradient.

qtot

total charge of the molecule.

dfac

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.

epssteep


if the norm of the gradient is greater than epssteep, a deepest-descent-step will be done.

epssearch


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

dqmax


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

mxls, dhls, ahls


parameters of linesearch:

ahls

start value

dhls

increment

mxls

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

transform


a switch for the transformation in the principal axis system.

lnumhess


a switch for the numerical Hessian.

lmd

a switch for a MD calculation.

Input Data Blocks Needed by UFF

$coord


cartesian coordinates of the atoms (default: $coord file=coord)

$ufftopology


contains a list of the next neigbours of each atom (see Section 21.2.4). Sometimes it is useful to enter the connectivity (in the input block nxtnei12 in the file ufftopology) by hand (not always necessary; default: $ufftopology file=ufftopology).

Beyond this uff reads the force field parameters for the atoms from the file parms.in. If this file exists in the directory from which one starts an uff calculation the program will use this file, if not the program reads the data from the file $TURBODIR/uff/parms.in. If one wants own atom types, one has to add these atoms types in the file parms.in. For each new atom type one has to specify the natural bond distance, the natural bond angle, the natural 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

$coord


contains the (updated) cartesian coordinates of the atoms (default: $coord file=coord).

$ufftopology


contains the full information of the topology of the molecule and the whole force field terms (see below; default: $ufftopology file=ufftopology).

$uffgradient


contains the accumulated cartesian analytical gradients (default: $uffgradient file=uffgradient).

$uffhessian


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.

21.2.5 Keywords for woelfling

Module WOELFLING reads options from data group

$woelfling

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

Furthermore,

 riter 0

counts the number of completed iteration (no option).

21.2.6 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 699 and for definition of the functionals the section 6.2 on page 240.

Example (default input):

$dft  
   functional b-p

gridsize integer or minteger


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

Example:

$dft  
   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 704). For the definition of gridtype 1 3, please refer to Eq. (16), (17) and, (19) in Ref. [201].

Example (default value):

$dft  
   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 [202], default is nkk 3. The usage of nkk makes sense only in the range 1 nkk 6.

Example:

$dft  
   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.

Example:

$dft  
  gridsize 9  
  ntheta 30  
  nphi 60

old_RbCs_xi


Original grids had not been carefully optimized for all atoms individually. This has now been done, which let to changes of ξ for Rb and Cs only resulting in minor improvements. If you have ongoing projects, which have been started with the old grids, you should continue using them with the keyword old_RbCs_xi.

Example:

$dft  
   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.  [201].

In addition, the reduction of spherical grid points near nuclei is suppressed, i.e. fullshell on is set (see page 673).

Note: the keyword radsize integer overrules the setting of incr; for more information, see p. 668.

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. 684 and p. 690).

Example (Recommendation):

$dft  
  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. [201]. With the following input,

$dft  
   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. [201]. Program stops after this testing.

reference


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:

$dft  
   functional b-p  
   reference

test-integ


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

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.

fullshell


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:

$dft  
   functional b-p  
   gridsize m4  
   fullshell

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 21.2.6 on page 684).

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 21.2.10 on page 733.

gridordering


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:

$dft  
   gridordering 0

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

Example:

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

$firstorder


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

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

$mom


This keyword enables the use of the maximum overlap method in unrestricted ridft calculations to access excited states self-consistently. [203]

$moprint


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 )

$point_charges


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:

thr=real


distance threshold for discarding redundant point charges, default value 10-6.

selfenergy


if given, the selfenergy of the point charge array will will be included in the energy and the gradient

nocheck


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

list


print all point charges in the output (default is to print the point charges only if less than 100 charges given)

$prediag


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

   dim(fock,dens)=6072  
   natoms=6  
   nshell=34  
   nbf(CAO)=108  
   nbf(AO)=98  
   dim(trafo[SAO<-->AO/CAO])=256  
   rhfshells=1

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

errvec=char

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


no DIIS

char= sFDs


use S-12FDS12 - transposed

Further suboptions:

maxiter=integer


maximum number of iterations used for extrapolation.

debug=integer


debug level (default: 0)

integer=1

print applied DIIS coefficients

integer=2

print DIIS matrix and eigenvalues, too

qscal=real


scaling factor in DIIS procedure: qscal > 1 implies some damping, qscal = 1.0: straight DIIS.

thrd=real


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.

Recommended:

errvec=sFDs leads to the following defaults:
qscal=1.2, for SCF runs: maxiter=6 and thrd=0.3, prediag is off; for DFT runs: maxiter=5 and thrd=0.1 prediag is on. If you want to switch off prediag put
$prediag none.

$scfdump


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

unit=integer


Fortran unit number for this file. Unit numbers 30,31, are recommended.

size=integer


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.

file=char


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:

$scfintunit  
 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

none


To perform a calculation without a start vector (i.e. use a core Hamiltonian guess).

file=char


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:

expanded


These MOs were obtained by projection form another basis set. They should not be used for wavefunction analysis.

scfconv=integer


The MOs are converged SCF MOs , the convergence criterion applied was 10-integer

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

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

noautomatic


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.

closedshell=real


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

individual


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

off

Do not perform integrals statistics

dscf

Perform integrals statistics for dscf

mpgrad

see mpgrad

dscf parallel 

see PARALLEL PROCESSING

Options dscf parallel, grad, mpgrad will be described in the related chapters.
If $statistics dscf has been given integral prescreening will be performed (which is an 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.

RHF/ROHF

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

$rohf


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.
$roothaan


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:

$roothaan  
    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 22.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.

UHF

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

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

$uhf


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

$uhfmo_beta


see $uhfmo_alpha

DFT
$dft  
 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:

b-lyp

b3-lyp

b3-lyp_Gaussian

(equivalent to the Gaussian98 keyword B3LYP with VWNIII)

bh-lyp

s-vwn

s-vwn_Gaussian

(equivalent to the Gaussian98 keyword SVWN with VWNIII)

tpss

tpssh

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.

RI

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.

$rij


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.

$ridft


Obsolete keyword - use $rij instead!

$rik


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

$ripop


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

RI-JK

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.

MARI-J

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.

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

The following options are available:

precision

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.

lmaxmom

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.

thrmom

Threshold for moment summation. For highly accurate calculations it should be set to 1 10-24.

nbinmax

number of bins per atom for partitioning of electron densities. Default value is 8 and hardly ever needs to be changed.

wsindex

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.

extmax

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:

$senex  
$dsenex  
 gridsize 5

LHF

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.

$dft  
 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)

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

asymptotic


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.

homo


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:

$soghf


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

$kramers


switches on Kramers-restricted formalism

$collinear


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

$gdiis


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.

$rx2c


switches on X2C calculation.

$rbss


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 [204].

$rlocal


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 [79].

All of these keywords are combinable with $soghf.

21.2.7 Keywords for Module riper

Many of the dscf and ridft keywords are also used by riper. To run the calculations, DFT method needs to be set in the $dft data group (see 21.2.6). The list of options specific to the riper module is presented below:

$riper  
  # CFMM control options  
  lmaxmom    20  
  nctrgt     10  
  wsicl      3.0  
  epsbext    1.0d-9  
  locmult    on  
  locmomadd  2  
  # PCG control options  
  lpcg       on  
  lcfmmpcg   on  
  lmxmompcg  20  
  pcgtol     1.0d-9  
  pcgtyp     sp  
  # remaining options  
  northol    5  
  lchgprj    on  
  kinescr    on  
  pqmatdiag  off  
  pqsingtol  1.0d-8

The following options are available:

#RI-CFMM options

lmaxmom

Maximum l-moment of multipole expansions for calculations of Coulomb interactions. Default value is 20 and hardly ever needs to be changed.

nctrgt

Target number of charges per lowest level box of an octree. Default value is 10. For very large systems its increase may lead to a moderate speed-up.

wsicl

The well-separateness criterion. The boxes with centers separated more than sum of their lengths times 0.5×wsicl are considered as well-separated. The standard value is 3.0 and hardly ever needs to be changed. The usage of wsicl makes sense only for values 2.0. For wsicl 3.0 increasing lmaxmom may be necessary.

epsbext

Specifies precision parameter for the multipole expansion, used to determine basis function extents. Default value is 1.010-9. Hardly ever needs to be changed.

locmult

If set to on, an additional local multipole scheme is used for the charges from non-well-separated boxes. 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 molecules and periodic systems, respectively.

locmomadd

If locmult is on increases the default order of local multipole expansions. Default value is 2 and probably never needs to be increased.

#PCG options

lpcg

If set to on the preconditioned conjugate gradient (PCG) algorithm is used for solving the RI equations. It is implemented for the molecular systems only. Default value is off.

lcfmmpcg

If lpcg is used, lcfmmpcg specifies whether the CFMM is applied for evaluation of the matrix-vector products needed for the iterative PCG solver. Not employing CFMM speeds up the calculations but halves the memory savings since full PQ matrix is stored. Default value is on.

lmxmompcg

Maximum l-moment of multipole expansions for calculations of Coulomb interactions within the PCG algorithm. Default value is 20. It should be set to the same or larger value than lmaxmom.

pcgtol

Sets the threshold parameter controlling accuracy of the PCG solver (see [99] 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 Chapter 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.

#remaining options

northol

Diagonalization of the Fock matrix and orthonormalization of orbital coefficients is done every northol SCF iteration. Default value is 5.

lchgprj

Specifies if charge projection is used for molecular systems. If set to off the output energy is identical as in case of ridft calculations for molecules.

kinescr

When set to on (default) the kinetic energy criterion is used for integral screening. Probably should never be changed.

pqmatdiag

If set to on full diagonalization of PQ matrix is performed when solving the RI equations. When diffuse basis functions are used the Cholesky decomposition may fail due to small eigenvalues. In this case the slower method based on a full diagonalization of the PQ matrix is necessary. Default value is off.

pqsingtol

If pqmatdiag is used pqsingtol sets threshold to remove small eigenvalues. Default value is 1.0 10-8.

If the $riper group is absent, the default values are used. For periodic systems, the following keywords are used to define the system:

$periodic n


Specifies the number of periodic direction. Value of n is 3 for a three-dimensional bulk, 2 for a two-dimensional surface slab and 1 for a one-dimensional system.

$cell


|a||b||c| α β γ

Unit cell parameters in form of six real numbers. 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. Lengths and angles in atomic units and degrees, respectively. Note that riper assumes that the cell vectors a and b are aligned along the x axis and on the xy plane, respectively. For 1D systems the periodic direction is defined by vector a. In case of 2D calculations first and second periodic directions are a and b, respectively.

$kpoints


nkpoints n1 n2 n3

Defines a Γ-centered mesh of points kα = (k1,k2,k3) with components

k = i∕n with i = - (n - 1)∕2,...,(n - 1)∕2
 j     j           j            j

in crystal coordinates.

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

periodic

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.

cell

Unit cell parameters in a form of six real values |a|, |b|, |c|, α, β, γ, where |a|, |b|, |c| are lengths of the appropriate cell vectors, α is the angle between vectors b and c, β is the angle between vectors a and c, and γ is the angle between vectors a and b. Default are atomic units and degrees. You can specify unit cell parameters in ┼ and degrees using cell ang.

content  
  label x y z  
end

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.

cluster  
  label x y z  
end

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.

charges  
  label charge  
end

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

ch_list  
  label charge  
end

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:

lmaxmom


Maximum order of the multipole expansions in periodic fast multipole method (PFMM). Default value is 25.

potval


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.

wsicl


Well-separateness criterion for PFMM. Default is 3.0.

epsilon


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

21.2.9 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 20. 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 20.

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:

$cosmo  
  epsilon=infinity  
  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  
  use_old_amat  
  use_contcav  
  no_oc

epsilon=real


defines a finite permittivity used for scaling of the screening charges.

allocate_nps=integer


skips the COSMO segment statistics run and allocates memory for the given number of segments.

no_oc


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)

use_old_amat


uses A matrix setup of TURBOMOLE 5.7

use_contcav


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:

$cosmo_atoms  
# 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.:

  SCREENING CHARGE:  
    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 20). If the $cosmo_isodens keyword is given without suboptions, a scaled iosodensity cavity with default settings will be created. Possible options are:

$cosmo_isodens


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

dx=real


spacing of the marching tetrahedron grid in ┼ (default: 0.3┼).

all_dens=real


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.

$cosmo_isorad


dx=real


spacing of the marching tetrahedron grid in ┼ (default: 0.3┼).

COSMO in MP2 Calculations: The iterative COSMO PTED scheme (see chapter 20) 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

$cosmo_correlated

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

Water:

h2o_25.pot

Ethanol:

ethanol_25.pot

Methanol:

methanol_25.pot

Tetrahydrofurane:

thf_25.pot

Acetone:

propanone_25.pot

Chloroform:

chcl3_25.pot

Tetrachloromethane:

ccl4_25.pot

Acetonitrile:

acetonitrile_25.pot

Nitromethane:

nitromethane_25.pot

Dimethylsulfoxide:

dimethylsulfoxide_25.pot

Diethylether:

diethylether_25.pot

Hexane:

hexane_25.pot

Cyclohexane:

cyclohexane_25.pot

Benzene:

benzene_25.pot

Toluene:

toluene_25.pot

Aniline:

aniline_25.pot

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

  SCREENING CHARGE:  
    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 20. 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.

21.2.10 Keywords for Modules grad and rdgrad

Many of the dscf and ridft keywords are also used by grad and rdgrad.

$drvopt

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.

21.2.11 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

$drvopt


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 RAM memory to be used by the run (here 50 MB), about 70% of available memory should be fine, because $maxcor specifies only the memory used to store derivatives of density and Fock matrices as well as the CPHF-RHS. Default is 200 MB.

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

$nosalc


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.

$noproj


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.

$nomw


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.

$isosub


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. 

$isosub  
 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:

$ironly


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

$ramanonly


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

$les


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. 

$les
 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.

$hotfcht


Triggers the generation of input files for hotFCHT (program to calculate Franck-Condon-factors by R. Berger and co-workers). See 14.4.

$sijuai_out


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.

21.2.12 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

21.2.13 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

or

$scfinstab rpat

or

$scfinstab urpa

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

$scfinstab ciss

or

$scfinstab cist

or

$scfinstab ucis

or

$scfinstab spinflip

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

$scfinstab soghf

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

$scfinstab singlet

or

$scfinstab triplet

or

$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  
  590  
  400 i

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

$soes  
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
$gw

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

$gw  
  nl 1  
  gam 0.001  
  output qpe.dat  
  rpa

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: qpe.dat. Output filename for the quasiparticle energies.

    rpa

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.

21.2.14 Keywords for Module rirpa

The keyword $rirpa allows to specify the following options,

npoints n


Number of frequency integration points, n (default is 60).

nohxx


HF energy calculation is skipped, (HXX = Hartree + eXact (Fock) eXchange).

rpaprof


Generates profiling output.

rpagrad


Switches on the gradients calculation for RI-RPA.

drimp2


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.

21.2.15 Keywords for Module egrad

egrad uses the same general keywords as escf and grad, see Sections 21.2.10 and 21.2.13.

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.

$nacme

flag to compute Cartesian non-adiabatic coupling vectors between the excited state of interest and the ground state [205]. This option requires the use of weight derivatives in section dft. It is only implemented for C1 symmetry.

21.2.16 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. If $maxcor is not found, its value is set to 200 MB.

$mp2energy


Calculation of MP2 gradient is omitted, only MP2 energy is calculated.

$freeze  
 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

$freeze  
 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 regarded not by mpgrad, moreover freezing of virtual orbitals is generally not supported by mpgrad.

All essential data groups for mpgrad may be generated by the preparation tool mp2prep, apart from $maxcor (see above) these are the following:

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

$mointunit  
 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.

$mp2pair


calculation of MP2 pair correlation energies.

21.2.17 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 21.2.6.

$cbas file=auxbasis


cross reference for the file specifying the auxiliary basis as referenced in $atoms. We strongly recommend using auxbasis sets optimized for the corresponding MO basis sets.

$freeze


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

$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 20


The data group $maxcor adjusts the maximum size of core memory in MB which will be allocated during the ricc2 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.

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

$laplace

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.

conv


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

$ricc2

ccs  
cis  
mp2     d1diag  
cis(d)  energy only  
cis(dinf)  
adc(2)  
cc2  
ccsd  
mp3  
mp4  
ccsd(t)  
restart  
norestart  
hard_restart  
nohard_restart  
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  
sos  
gsonly  
d1diag  
intcorr

specifies the ab initio models (methods) for ground and excited states and the most important parameters and thresholds for the solution of the cluster equations, linear response equations or eigenvalue problems. If more than one model is given, the corresponding calculations are performed successively. Note: The CCS ground state energy is identical with the SCF reference energy, CCS excitation energies are identical to CIS excitation energies. The MP2 results is equivalent to the result from the rimp2 module. cis(dinf) denotes the iterative CIS(D) variant CIS(D). The option ccsd(t) request a CCSD calculation with the perturbative triples correction, CCSD(T), and as a side result also the CCSD[T] energy will be printed.

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.

(no)restart


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

(no)hard_restart


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.

conv

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.

oconv


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.

lindep


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.

maxiter


gives the maximum number of iterations for the solution of the cluster equations, eigenvalue problems or response equations (default: 25).

mxdiis


is the maximum number of vectors used in the DIIS procedures for ground state or excitation energies (default: 10).

maxred


the maximum dimension of the reduced space in the solution of linear equations (default: 100).

iprint


print level, by default set to 1 or (if given) the the value of the $printlevel data group.

fmtprop


Fortran print format used to print several results (in particular one-electron properties and transition moments) to standard output.

geoopt


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.

scs


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.

sos


the SOS parameters cos=1.3 and css=0.0 are applied. This keyword can presently only be used in connection with MP2.

d1diag


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

intcorr


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.

$rir12

ansatz  
r12model  
comaprox  
cabs  
examp  
r12orb  
pairenergy  
corrfac  
cabsingles  
f12metric

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.

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.

no_f12metric

f12metric


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.

$excitations

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  
xgrad    states=(ag{3} 1)  
conv    = 6  
thrdiis = 2  
preopt  = 3  
leftopt  
bothsides  
oldnorm

In this data group you have to give additional input for calculations on excited states:

irrep


the irreducible representation.

multiplicity


spin multiplicity (1 for singlet, 3 for triplet); default: singlet, not needed for UHF.

nexc

the number of excited states to be calculated within this irrep and for this multiplicity.

npre

the number of roots used in preoptimization steps (default: npre = nexc).

nstart


the number of start vectors generated or read from file (default: nstart = npre).

spectrum


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.

tmexc


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.

exprop


require calculation of first-order properties for excited states. For the states option see spectrum option above; for details for the operators input see below.

xgrad


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.

conv

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.

preopt


convergence threshold used for preoptimization of CC2 eigenvectors is set to 10-preopt (default: 3).

thrdiis


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

leftopt


If the flag leftopt is set the left eigenvectors are computed (default is to compute the right eigenvectors, for test purposes only).

bothsides


The bothsides flag enforces the calculation of both, the left and the right eigenvectors (for test purposes only).

bothsides


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.

$response

fop unrelaxed_only operators=diplen  
sop operators=(diplen,diplen) freq=0.077d0  
gradient  
conv = 6  
zconv = 6  
semicano  
nosemicano  
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:

fop

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

sop

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.

gradient


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.

conv

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.

zconv


convergence threshold for the norm of the residual vector in the solution of the Z vector equations will be set to 10-zconv.

semicano


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

nosemicano


use non-canonical formulation for the calculation of (transition) one-electron densities. Default is to use the semi-canonical formulation.

thrsemi


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

zpreopt


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.

nozpreopt


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:

operators=diplen,dipvel


input of operator labels for first-order properties, transition moments, etc. Currently implemented operators/labels are

overlap

overlap (charge) operator: the integrals evaluated in the AO basis are μ|ν

diplen

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

dipvel

dipole operator in velocity gauge: μ|∇i|ν
(all three components, individual components can be specified with the labels xdipvel, ydipvel, zdipvel).

qudlen

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⟩
     2

(including nuclear contributions) is given.

angmom

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

nef

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

states=all


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.

$D2-diagnostic


Calculate the double-substitution-based diagnostics D2.

$cc2_natocc


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

21.2.18 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 21.2.6.

$cbas file=auxbasis


Auxiliary basis set for RI approximation. For details Section 21.2.17.

$freeze


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

$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 5000


The data group $maxcor adjusts the maximum size of core memory in MB 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.

$laplace

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.

$pnoccsd

mp2  
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  
sos

mp2

specifies the ab initio model (method). The current release version is restriced to MP2 which is also the default model.

localize

specifies the localization method; possible choices are boys for Foster-Boys, pm for Pipek-Mezey, and none for canonical (deprecated, only meant for testing) orbitals. Default are Foster-Boys orbitals. (Pipek-Mezey orbitals become very delocal with diffuse basis sets.)

prepno

switch between to algorithms for the generation of OSV coefficients. Possible choices are full for an O(N4) scaling direct diagonalization (cmp. [206]) and davidson for an O(N3) scaling iterative diagonalization (cmp.  [20]). davidson is the default choice and recommended.

mxrdim

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. [206]). Usually there is no need to touch this parameter.

tolpno

specifies the PNO truncation threshold. Default value: 10-7.

tolosv

specifies the OSV truncation threshold. If no given tolosv is set to 0.1×tolpno, which is the recommended value.

tolri

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.

tolpair

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

opnos

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

tolcapno

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.

tolosc

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.

tolopno

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.

toloso

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.

projector

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.

conv

The conv parameter gives the convergence threshold for the ground state energy as 10-conv. The default threshold is 10-7.

oconv

The oconv parameter gives an additional threshold for the residual of the ground state equations as < 10-oconv. The default threshold is 10-3.

lindep

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.

maxiter

gives the maximum number of iterations for the solution of the cluster equations, eigenvalue problems or response equations (default: 25).

mxdiis

is the maximum number of vectors used in the DIIS procedures for the ground state equations (default: 10).

maxred

not used in the current release.

scs

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.

sos

the SOS parameters cos=1.3 and css=0.0 are applied.

$rir12

for the description of this data group see Sec. 21.2.17.

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

logarithm


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

scale


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

NOTES:
  • 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:

maxiter=n


maximum number of iterations for the iterative conversion procedure internal cartesian coordinates (default: 25).

qconv


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:

none

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

numgeo=integer

number of structures used

maxgeo=integer

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

ingeo=integer

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

NOTES: 

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:

numpul=integer


number of geometries to be utilized

maxpul=integer


maximum number of geometries

minpul=integer


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

fail=real


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

diagonal

 
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.

offreset

 
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!

allow=real

 
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)

scale=real

 
scaling factor for the input hessian (default: 1.0).

threig=real

 
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:

reseig=real

default: texttt0.005.

thrbig=real

 
upper bound for eigenvalues of the hessian; if any eigenvalue exceeds thrbig, it will limited to this value (default: 1000.0).

damping=real

 
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:

on/off


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!

diag=suboptions


provide a diagonal force constant matrix with:

available suboptions are:

real


this will lead to an assignment of diagonal elements (default: 1.0)).

default

 
this will lead to an assignment of initial force constant diagonals depending on the coordinate type.

individual

 
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

carthess

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

$grad

 
Cartesian coordinates and gradients calculated in subsequent optimization cycles. Entries are accumulated by one of the gradient programs (grad, mpgrad, rimp2, ricc2, egrad, etc.).

$egrad

 
Basis set exponents scale factors and their gradients as calculated in subsequent optimization cycles. Entries are accumulated by one of the gradient programs.

$globgrad

 
Global scale factors and gradients as calculated in subsequent optimization cycles. Entries are accumulated by the grad or aoforce program.

$corrgrad

 
Allows to augment internal SCF gradients by approximate increments obtained from treatments (e.g. correlation or relativistic) on higher level. See the example below.

$corrgrad  
#  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:

format=format

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

$coord


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.

$basis


updated basis set exponents, basis sets contraction coefficients or scaling factors, if $optimize basis on has been specified.

$global


updated global scaling factor for all basis set exponents, if $optimize global on has been specified.

$forceapprox


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

$forcestatic


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

$oldgrad


Then the relax program accumulates all subsequent coordinates and gradient as used in optimization in this data group (or a referenced file). This overrides the input of old coordinate and gradient data from data blocks $grad, $egrad, …as accumulated by the grad program.

degrees

21.2.20 Keywords for Module statpt

$statpt  
  itrvec      0  
  update      bfgs  
  hssfreq     0  
  keeptmode  
  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:

$statpt  
  itrvec 0

Following options are available:

itrvec


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.

update


Method of hessian update. For minimization default is BFGS, for TS search default is Powell and none is for no update.

hssfreq


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!

keeptmode


Freezing transition vector index.

hssidiag


diagonal hessian elements for diagonal Hessian guess (default: 0.5).

radmax


Maximum allowed value for trust radius (default: 0.3).

radmin


Minimum allowed value for trust radius (default: 1.0d-4).

tradius


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

Convergence criteria  

threchange

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

thrmaxdispl

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

thrmaxgrad

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

thrrmsdispl

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

thrrmsgrad

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

All values are in atomic units.

21.2.21 Keywords for Module moloch

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

$properties  
    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.

moments


if moment is active you need

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

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

potential


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.

localization


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

$mulliken  
 spdf molap netto irpspd irpmol mommul

to perform a Mulliken population analysis. The options specify the output data:

spdf

print molecular orbital contributions to atomic s, p, d,…-populations

molap

print molecular orbital contributions to overlap populations

netto

print atomic netto populations

irpspd

print contributions of (irreducible) representations to atomic s,p,d,…-populations

irpmol

print contributions of (irreducible) representations to overlap populations

or

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

or

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

momao

print MO contributions to occupation numbers of modified atomic orbitals (MAOs).

maodump

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

symb

atom symbol

list

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

nmao=i


means number of MAOs to be included

method=meth


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.

threshold=r


means the threshold to be applied for the selection criteria occ or eig (default: 0.1).

Example:

$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


plot


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.

fit


if fit is active you need

$vdw_fit  
shell     number_of_gridpoints     distance_from_vdW_surface  
refine    value_of_potential

shell

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

refine

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.

$vdw_radii  
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].

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

$mvd


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

$moments


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.

$pop


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:

lall

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.

thrpl=real


Contributions smaller than thrpl are not displayed (default: 0.01).

overlap

Mulliken atomic overlap matrix is displayed.

netto

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

Example:

$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 [163]. 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.

Example:

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

$nbonmb  
  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 [164] 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:

eig


MAOs are sorted according to their eigenvalue (those with largest EW finally are chosen). This is the default.

occ


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:

fix


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.

max


Maximum of numbers calculated from fix and thr=0.1 is taken.

mix


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:

info


Eigenvalues and occupations for each MAO are written to output.

dump


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.

$localize


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

sweeps=integer


maximum number of orbital rotations to get LMOs; default value is 10000 (sometimes not enough, in particular for highly delocalised systems).

thrcont=real


lower threshold for displaying MO and Mulliken contributions (default: 0.1).

CAO


LMOs are written to file in the CAO basis (instead of AO).

pipmez


Pipek-Mezey localization is used.

$wfn


triggers the generations of a wfn file. It can be used in dscf/ridft single-point calculations or in ricc2/egrad gradient calculations.

$esp_fit


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:

$esp_fit
shell i1 s1
shell i2 s2
...

Integer numbers i define the number of points for the respective shell, real numbers s constants added to radii (default corresponds to one shell with s=1.0).

A parameterization very close to that by Kollman (U.C. Singh, P.A. Kollman, J. Comput. Chem. 5(2), 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.

$pointval


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

pot

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.

fld

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

dens

has to be set, if additionally to one of the above quantities also the density is to be computed.

xc

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:

xyz

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.

plt

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.

map

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.

txt

a format compatible with gOpenMol for visualization of vectors v. The format is x,y,z,vx,vy,vz.

vec

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.

cub

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:

$pointval  
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:

$pointval  
  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)
(21.1)

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

$md_status


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:

$md_status  
  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

$title


Arbitrary title

$log_history  
 100                mdlog.P  
 71                 mdlog.Q

$ke_control  
  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.

type


can be one of orth, elps, or none, for orthorhombic, ellipsoidal, or no barrier (the default) respectively.

limits


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.

constant


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.

springlen


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

temperature


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

$constraints


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.

tolerance


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.

adjpercyc


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:

$md_action  
     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:

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

$md_action  
     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.

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

$md_action  
     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.

$surface_hopping

This keyword allows to carry out Tully-type fewest switches surface hopping (SH) [207]. This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [208]; for the current implementation see Tapavicza et al. 2011 [209]. 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

$md_action  
 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:

$nacme


needed to compute non-adiabatic couplings; this keyword requires the use of weight derivatives in section dft

$nac


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 [210,211]. 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.

21.2.24 Keywords for Module mpshift

In order to control the program execution, you can use the following keywords within the control file:

$csmp2


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.

$mointunit


Scratch file settings for an MP2 calculation. Please refer to Section 21.2.16 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 21.2.6);
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.

21.2.25 Keywords for Parallel Runs

On all systems the parallel input preparation is done automatically. Details for the parallel installation are given in Section 3.2.1. The following keywords are necessary for all parallel runs:

$parallel_platform architecture

Currently the following parallel platforms are supported:

SMP

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

MPP

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

cluster

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

SGI

similar to SMP, but here the server task is treated differently: the MPI implementation on the SGIs would cause this task to request too much CPU time otherwise.

If you want to run mpgrad, $traloop has to be equal to or a multiple of the number of parallel workers.

For very large parallel runs it may be impossible to allocate the scratch files in the working directory. In this case the $scratch files option can be specified; an example for a dscf run is given below. The scratch directory must be accessible from all nodes.

$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 all programs employing density functional theory (DFT) (i.e. dscf/gradand ridft/rdgrad) $pardft can be specified:

$pardft  
       tasksize=1000  
       memdiv=0

The tasksize is the approximate number of points in one DFT task (default: 1000) and memdiv says whether the nodes are dedicated exclusively to your job (memdiv=1) or not (default: memdiv=0).

For dscf and grad runs you need a parallel statistics file which has to be generated in advance. The filename is specified with
$2e-ints_shell_statistics    file=DSCF-par-stat
or
$2e-ints’_shell_statistics    file=GRAD-par-stat
respectively.

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

$statistics  dscf parallel  
$2e-ints_shell_statistics    file=DSCF-par-stat  
$parallel_parameters  
       maxtask=400  
       maxdisk=0  
       dynamic_fraction=0.300000

and for a grad statistics run:

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

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.

For parallel grad and rdgrad runs one can also specify:

$grad_send_dens

This means that the density matrix is computed by one node and distributed to the other nodes rather than computed by every slave.

In the parallel version of ridft, the first client reads in the keyword $ricore from the control file and uses the given memory for the additional RI matrices and for RI-integral storage. All other clients use the same amount of memory as the first client does, although they do not need to store any of those matrices. This leads to a better usage of the available memory per node. But in the case of a big number of auxiliary basis functions, the RI matrices may become bigger than the specified $ricore and all clients will use as much memory as those matrices would allocate even if that amount is much larger than the given memory. To omit this behavior one can use:

$ricore_slave integer

specifying the number of MBs that shall be used on each client.

For parallel jobex runs one has to specify all the parallel keywords needed for the different parts of the geometry optimization, i.e. those for dscf and grad, or those for ridft and rdgrad, or those for dscf and mpgrad.