TURBOMOLE input is keyword-directed. Keywords start with a ’$’, e.g. $title. Comments may be given after $dummy, or by a line starting with #; these lines are ignored by TURBOMOLE. Blank lines are also ignored. Keywords may be in any order unless stated otherwise below.
The sample inputs given below should help to give an idea how the keywords are to be used. Complete control files are provided in Chapter 21. An alphabetical list of all keywords is given in the index.
The control file and other input files references therein must end with this keyword:
General information defining the molecular system: nuclear coordinates, symmetry, basis functions, number of occupied MOs, etc. which are required by every module.
give title of run or project here.
Schönflies symbol of the point group. All point groups are supported with the
exception of NMR shielding and force constant calculations etc. which do not work
for groups with complex irreps (C3, C3h, T, etc). Use a lower symmetry group in this
case.
Example:
the basis set
and the auxiliary (fitting) basis for RIDFT calculations
the ECP if this is used.
The files basis, ecp and jbas must provide the necessary information under the labels specified in $atoms.
This data group specifies the number of Cartesian components of basis functions (i.e. 5d
and 7f in 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).
Specification of MO occupation for RHF, e.g.
MO occupation of open shells and number of open shells. type=1 here means that there is
only a single open shell consisting e.g. of two MOs:
Roothaan parameters for the open shell, here a triplet case. define recognizes most cases
and suggests good Roothaan parameters.
For further information on ROHF calculations, see the sample input in Section 21.6 and the tables of Roothaan parameters in Section 6.3.
directs the program to carry out a UHF run,e.g.
The specification of MO occupation for UHF, $uhf overwrites closed-shell occupation specification.
Most post-SCF programs (aoforce, ccsdf12, egrad, escf, evib, mpgrad, mpshift, pnoccsd, ricc2, rirpa) in TURBOMOLE read the data group
to control the amount of dynamically allocated memory. The integer number (above “500”) one can (optionally) give a SI prefix (kb, mb, gb) or an IEC prefix (kib, mib, gib) to specify the unit. If no unit is specified mib (binary megabytes, i.e. 10242 bytes) are used.
In addition one can specify a reference for parallel calculations:
If not given the specified memory will be used per thread. Note that in release versions before V7.2 the memory was specified per process. For sequential calculations this sub-option has no effect.
If $maxcor is not given 500 MiB per thread will be used as limit for dynamically allocated memory.
Note:
The four keywords above are set by define, but are not necessary.
$statistics dscf
or
$statistics mpgrad
Only a statistics run will be performed to determine file space requirements as specified for dscf or mpgrad. On return the statistics option will be changed to $statistics off.
means current step. Keyword and data group (as e.g. dscf) is set by every program and removed on successful completion.
Keyword and data group (as e.g. relax) set by every program on successful completion.
General file cross-references:
It is convenient not to include all input in the control file directly and to refer instead to other files providing the corresponding information. The above cross references are default settings from define; you may use other file names. define will create most of these files. Examples of these files are given below in the samples.
contains atom specification—type and location—and the bonds and internal
coordinates convenient for geometry optimizations.
specification of basis sets.
specification of effective core potentials.
MO vectors of SCF or DFT calculations for RHF or UHF runs.
keywords and data groups set by unrestricted dscf or ridft runs. Contain natural
MO vector and orbital occupation.
energies and gradients of all runs, e.g. for documentation in a geometry optimizations.
approximate force constant for geometry optimizations.
With the parameters in $redund_inp the generation of redundant internal coordinates can be modified. All entries have to be made in the control file before invoking the ired option. Important options are:
print parameter for debug output: The larger n is, the more output is printed
n ≥ 0,n ≤ 5 (default: 0)
method for generating and processing of redundant internal coordinates
n ≥ - 3,n ≤ 3,n ≠ 0 (default: 3)
Values for the metric option:
“Delocalized Coordinates”
The BmBt matrix is diagonalized for the complete set of redundant internal
coordinates, matrix m is a unit matrix.
Delocalized Coordinates obtained with a modified matrix m, the values of m can be defined by user input (see below).
“Hybrid Coordinates”
Natural internal coordinates are defined as in the old iaut option. If a cage
remains, delocalized coordinates (as for n=1) are defined for the cage.
Very simular to the n = 1 option, but for the remaining cage delocalized coordinates with modified matrix m are defined as for n = -3.
“Decoupled coordinates”
The redundant coordinates are divided into a sequence of blocks. These
are expected to have decreasing average force constants, i.e. stretches,
angle coordinates, torsions and “weak” coordinates. The 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.
“Generalized natural coordinates”
Natural internal coordinates are defined first, for the remaining cage
decoupled coordinates are defined.
a positive real number, which is an approximate “force constant”, can be read in for each
type of coordinate (see below). The force constants are used for the definition of the matrix
m in BmBt.
The matrix m is assumed to be a diagonal matrix. For each type of coordinate a different value for the force constants mii can be read in. Types of coordinates are:
bond stretch (default: 0.5)
inverse bond stretch (default: 0.5)
bond angle (default: 0.2)
Out of plane angle (default: 0.2)
dihedral or “torsional” angle (default: 0.2)
Special angle coordinate for collinear chains, bending of the chain a–b–c in the plane of b–c–d (default: 0.2)
bending of the chain a–b–c perpendicular to the plane of b–c–d
(default: 0.2)
stretch of a “weak” bond, i.e. the bond is assumed to have a very low force constant,
e.g. a “hydrogen bond” or a “van der Waals bond”
(default: 0.05)
inverse stretch of a weak bond (default: 0.05)
bond angle involving at least one weak bond (default: 0.02)
Out of plane angle for weak bonds (default: 0.02)
dihedral angle for weak bonds (default: 0.02)
linc coordinate for weak bonds (default: 0.02)
linp coordinate for weak bonds (default: 0.02)
One has to specify only the Cartesian coordinates (data group $coord) to start a uff run. The program uff requires the data groups $uff, $ufftopology, $uffgradient and $uffhessian. If these keywords do not exist in the control file the program will generate these data groups.
The data group $uff contains the parameters described below. The default values—in the control file—are:
The explanation of the variables are as follows:
number of max. optimization cycles (maxcycle=1: single-point calculation).
can have the values +1 or -1. If modus = -1 only the topology will be calculated.
each nqeq cycle the partial charges will be calculated. If nqeq = 0, then the partial charges are calculated only in the first cycle, if the file ufftopology does not exist.
switch for the different types of force field terms:
bond terms will be calculated.
angle terms will be calculated.
torsion terms will be calculated.
inversion terms will be calculated.
non bonded van der Waals terms will be calculated.
non bonded electrostatic terms will be calculated.
convergence criteria for energy and gradient.
total charge of the molecule.
distance parameter to calculate the topology. If the distance between the atoms I and J is less than the sum of the covalent radii of the the atoms multiplied with dfac, then there is a bond between I and J.
if the norm of the gradient is greater than epssteep, a deepest-descent-step will be
done.
if the norm of the gradient is smaller than epssearch, no line-search step will be done after
the Newton step.
max. displacement in a.u. for a coordinate in a relax step.
parameters of linesearch:
start value
increment
number of energy calculations
modification parameter for the eigenvalues of the Hessian (see below):
f(x) = x * (alpha + beta* exp(-gamma* x)).
a switch for the transformation in the principal axis system.
a switch for the numerical Hessian.
a switch for a MD calculation.
cartesian coordinates of the atoms (default: $coord file=coord)
contains a list of the next neigbours of each atom (see Section 20.2.5). Sometimes it is
useful to enter the connectivity (in the input block nxtnei12 in the file ufftopology)
by hand (not always necessary; default: $ufftopology file=ufftopology).
Beyond this uff reads the force field parameters for the atoms from the file parms.in. If this file exists in the directory from which one starts an uff calculation the program will use this file, if not the program reads the data from the file $TURBODIR/uff/parms.in. If one wants own atom types, one has to add these atoms types in the file parms.in. For each new atom type one has to specify the natural bond distance, the natural bond angle, the natural 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.
contains the (updated) cartesian coordinates of the atoms (default: $coord
file=coord).
contains the full information of the topology of the molecule and the whole force field
terms (see below; default: $ufftopology file=ufftopology).
contains the accumulated cartesian analytical gradients (default: $uffgradient
file=uffgradient).
contains the cartesian analytical Hessian;
(default: $uffhessian file=uffhessian0-0).
The 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:
|
Here are I and J the number of the atoms, d the distance in a.u. and BO is the bond order.
The angle terms follow, starting with the number of the angle terms. In each line is one angle term:
|
Here are J,I and K the atoms number, where atom I is in the apex. “wtyp” is the angle type and has the following values:
linear case
trigonal planar case
quadratic planar case
octahedral case
all other cases.
θ is the angle value in degree. 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:
|
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:
J (sp3)–K (sp3)
like ttyp=1, but one or both atoms are in Group 16
J (sp2)–K (sp3) or vice versa
like ttyp=2, but one or both atoms are in Group 16
like ttyp=2, but J or K is next a sp2 atom
J (sp2)–K (sp2)
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 and L are the atom numbers. Atom I is the central one. ityp1, ityp2, ityp3 are the types of the inversions:
atom I is C and atom L is O
like ityp=10, but L is any atom
I is P
I is As
I is Sb
I is Bi
all other cases.
ω1,ω2 and ω3 are the values of the inversion angles in degree.
The nonbond terms follow starting with the number of the non-bonded terms. In each line is one nonbond term:
|
Here I and J are the atom numbers, d the distance in a.u. Then the partial charges follow.
If the determination of the molecule connectivity failed, you can specify the block nxtnei12 in the file ufftopology. Then the program calculates the other blocks.
Based on the numbers of the next neighbours (block nxtnei12 in the file ufftopology) the program tries to determine the UFF type of an atom. The following rules are implemented: If the atom has three next neighbours and it is in the nitrogen group, then it has a hybridization three. If it is not in the nitrogen group, it has hybridization two. If the atom has four next neighbours and it is in the carbon group, it has hybridization three. If it is not in the carbon group, it becomes hybridization four. If the number of next neighbours is six, then it gets the hybridization six.
Since the smallest eigenvalues λi of the Hessian has the greatest influence on the convergence of the geometry optimization, one can shift these values with
|
and calculates a new Hessian with these modified eigenvalues.
Module WOELFLING reads options from data group
The below values of the options are default values with the following meaning:
Number of interpolated structures for optimization.
Number of input structures provided by user.
Align input structures by translation/rotation 0=yes, 1=no.
Maximum number of iterations.
Threshold for accuracy of LST-interpolation.
Threshold for mean of norms of projected gradients.
Use standard optimization from initial LST-path (method q) or grow reaction path (method qg).
Furthermore,
counts the number of completed iteration (no option).
Convergency criterion for the root mean square of the density matrix. If you want to
calculate an analytical MP2 gradient (program mpgrad) real should be 1.d-7 or less.
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.
Specification of the functional, default is BP86, printed as functional b-p. For
all possible—and useful—functionals, please refer to page 765 and for definition
of the functionals the section 6.2 on page 250.
Example (default input):
Specification of the spherical grid (see section 20.2.7 on page 765). Default is
gridsize m3.
Example:
—not recommended for use—
Specification of the mapping of the radial grid.
Possible values for gridtype are 1,…,6, but gridtype 4 to 6 is only for the use with functional lhf (see page 770). For the definition of gridtype 1 – 3, please refer to Eq. (16), (17) and, (19) in Ref. [232].
Example (default value):
Flag for debugging. debug 0 means no debug output (default), debug 1 means some
output, debug 2 means a lot more output. Be careful!
Specification of the sharpness of the partition function as proposed by Becke [233],
default is nkk 3. The usage of nkk makes sense only in the range 1 ≤ nkk
≤ 6.
Example:
—not recommended for use—
Only for 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:
Original grids had not been carefully optimized for all atoms individually. This has
now been done, which let to changes of ξ for Rb and Cs only resulting
in minor improvements. If you have ongoing projects, which have been
started with the old grids, you should continue using them with the keyword
old_RbCs_xi.
Example:
Specifies the number of radial grid points. Default values depend on type of atom and
grid (see keyword gridsize). The formula for the radial gridsize is given
as,
|
ioffrad is 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).
Serves to increase quadrature grids; this is recommended in case of very
diffuse wavefunctions. With the keyword diffuse grids are modified by
changing the linear scaling factor ξ of radial grid points and the radial
gridsize:
radsize radsize + incr
ξ ξ * scal
diffuse integer | 1 | 2 | 3 | 4 | 5 | 6 |
incr | 1 | 2 | 3 | 4 | 5 | 6 |
scal | 1.5 | 2.0 | 2.8 | 4.0 | 5.0 | 6.0 |
For information about the linear scaling parameter ξ, see Eq. (16)–(19) and Table 1 in Ref. [232].
In addition, the reduction of spherical grid points near nuclei is suppressed, i.e. fullshell on is set (see page 741).
Note: the keyword radsize integer overrules the setting of incr; for more information, see p. 736.
Recommendation: For diffuse cases use gridsize m4 (or larger) in combination with diffuse 2 and check the number of electrons; for more difficult cases use diffuse 4. In case of doubt, verify the calculation with a larger grid, i.e. gridsize 7.
The test suite example $TURBODIR/TURBOTEST/dscf/short/H3PO4.DSCF.DIFFUSE provides an example of usage; this also gives reasonable values for damping and orbitalshift to reach convergence in this and similar cases, see $scfdamp and $scforbitalshift (p. 750 and p. 756).
Example (Recommendation):
—for developers only—
Radial grid points have a linear scaling parameter ξ, see Eq.(16)–(19) and Table 1 in
Ref. [232]. With the following input,
one performs a numerical integration for the density and the exchange correlation term for ξ = 0.5,(0.01),2.0 for given MOs and functional. NOTE: only molecules with a single atom type can be used. The results serve to establish stable, optimal ξ values, see Figure 1 in Ref. [232]. Program stops after this testing.
Usage of the reference grid, which is a very fine grid with very tight thresholds. The
default values for the different variables are:
Please refer to the corresponding sub-keywords for explanation.
If you want to use the reference grid, you have to skip the keyword gridsize, and type instead reference. Example:
Check if the selected grid is accurate enough for the employed basis-set by performing
a numerical integration of the norm of all (occupied and virtual) orbitals. Useful for
LHF. 771.
Grid points are sorted into batches, which are then processed. This increases
efficency. This should be changed only by developers. Default is batchsize
100.
Standard grids have reduced number of spherical grid points near nuclei. With the
keyword fullshell this reduction is suppressed. Reference grid (see keyword
reference) always has full spherical grids with 1202 points. Should be used to
checked the influence of spherical grid reduction.
Example for the usage of fullshell:
—for developers only—
Values of real effects efficiency of the quadrature, default is symblock1 0.001 and
symblock2 0.001, one can try higher or smaller values.
—not recommended for use—
Where xparameter (default) can be: sgrenze (8), fgrenze (10), qgrenze (12),
dgrenze (12) and, fcut (14). These parameters control neglect of near zeros of
various quantities. With xparameter integer one changes the default. integer
larger than defaults will increase the numerical accuracy. Tighter threshold
are set automatically with keyword $scfconv (see section 20.2.7 on page
750).
Includes the derivatives of quadrature weights to get more accurate results. Default is
that the derivatives of quadrature weights will be not considered, see section 20.2.12
on page 811.
Grid points are ordered into batches of neighbouring points. This increases efficiency,
since now zeros can be skipped for entire batches. gridordering is default for serial
version, not for the parallel one. You cannot use weight derivatives and
gridordering together.
Example for switching off gridordering:
Within the semi-numerical integration scheme employed for local hybrid functionals,
S- and P-junctions feature pre-screenings of the analytical two-center integrals with
respect to shell pair overlap and shell pair conjunction through the density matrix,
respectively. This allows to neglect certain shell pairs and thus accelerates the
calculation. s-junc and p-junc set the corresponding thresholds to 10-integer.
Smaller values increase the number of neglected shell pairs and thus reduce the
calculation cost, but also decrease the accuracy of the results. Note that for gradient
and scf calculations the P-junctions also depend on the grid weights. If
larger grids are employed seeking more accurate results, the threshold for
P-junctions should thus be adjusted (increasing the value to p-junc 6). The
defaults settings are s-junc 6 and p-junc 6 for grad, rdgrad, and ridft.
Too small values for s-junc and p-junc might result in non-converging
calculations or deteriorated results. This holds particularly for excited state
calculations. For escf the default values are therefore s-junc 8 and p-junc
8.
Specification of external electrostatic field(s). The specification may take place either by
Ex, Ey, Ez or by x, y, z, |E|. See also $fldopt.
Example:
Requests calculation of occupation numbers at a finite temperature T. For an orbital with
the energy εi the occupation number ni ∈ is calculated as
where μ is the Fermi level. The factor f = 4k∕ is chosen to yield the same slope at the Fermi level as the Fermi distribution.
Calculation of the fractional occupation numbers starts when the current 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: 1⋅10-3). This prevents oscillations at the end of the SCF procedure.
Calculation of fractional occupation numbers often requires much higher damping and orbital shifting. Please adjust the values for $scfdamp and $scforbitalshift if you encounter convergence problems.
In UHF runs this option can be used to automatically locate the lowest spin state. In order to obtain integer occupation numbers tmend should be set to relatively low value, e.g. 50 K.
Calculation of fractional occupation numbers should be used only for single point calculations. When used during structure optimizations it may lead to energy oscillations.
The optional nue value (number of unpaired electrons) can be used to force a certain multiplicity in case of an unrestricted calculation. nue=0 is for singlet, nue=1 for dublet, etc.
The option addTS adds the entropic contribution of the smearing to the total energy which is important for very high temperatures only.
Finally the option noerf will use the full correct Fermi statistics rather than the term given above.
Perform first-order SCF-calculation, i.e. perform only one SCF-iteration with the start MOs
(which should be the orthogonalized MOs of two independent subsystems as is explained in
detail in Chapter 18).
Specification of options related with external electrostatic fields. The following options are
available:
Calculate numerical 1st derivative of SCF energy with respect to electrostatic
field (default: off), increment for numerical differentiation is edelt (see below).
Calculate numerical 2nd derivative of SCF energy with respect to electrostatic
field (default: off), increment for numerical differentiation is edelt.
Increment for numerical differentiation (default: 0.005).
Calculate SCF energy for non-zero external electrostatic
fields defined in $electrostatic field.
Calculate SCF energy for one external field definition and dump disturbed
MOs onto $scfmo. This enables to evaluate properties or perform geometry
optimizations in the presence of an external field.
Caution: don’t use the RI approximation for all these calculations since this will lead to non-negligible errors!!
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:
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.
If this keyword is set the energies and symmetry labels of all occupied MOs will be dumped
to this data group. This may be helpful to draw mo-diagrams. If only has been set only the
start MOs are dumped and the program quits.
nirreps will hold the total number of displayed orbitals after the successful run.
This keyword enables the use of the maximum overlap method in unrestricted ridft
calculations to access excited states self-consistently. [234]
If this keyword is present all occupied orbitals are dumped to standard output. Be
careful about this option as it can create huge output files in case of many basis
functions.
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))
This data group will be written after an UHF calculation (together with $natural orbital
occupation) and contains the natural space orbitals (same syntax as $scfmo).
This data group will be written after an UHF calculation (together with $natural
orbitals) and contains the occupation of natural orbitals (syntax as any data group
related with orbital occupation information, e.g. $closed shells), e.g.
concerns the first SCF iteration cycle if start MOs from an EHT guess are used.
The SCF iteration procedure requires control mechanisms to ensure (fast) convergence, in TURBOMOLE these are based on orbital energies ϵi of the preceeding iteration used for level shifting and damping (besides DIIS, see below). This feature cannot be used in the first iteration if EHT MOs are employed as start, since ϵi are not available. The keyword $prediag provides ’ϵi of the zeroth iteration’ by diagonalization of occ–occ and virt–virt part of the first Fock matrix, to allow for level shifting etc.. See $scfdiis below.
Try a dscf restart. The program will read the data group $restartd (which must exist,
also $scfmo has to exist!) and continue the calculation at the point where it ended before. If
the additional option twoint is appended, the program will read the 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.
Data provided by a previous dscf run that has been interrupted. This keyword is created
when $scfdump was given.
is set by define so usually no changes are necessary. The dimensions must be greater or
equal to those actually required, i.e. you can delete basis functions and keep rundimensions.
This keyword is not necessary for small cases.
Example:
SCF convergency criterion will be 10-integer for the energy. Gradients will only be evaluated
if integer > 6.
Damping parameters for SCF iterations in order to reduce oscillations. The old
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.
Flags for debugging purposes. Following options are available:
Output level concerning molecular orbitals. integer=0 (default) means minimal
output, >1 will output all start MOs and all MOs in each iteration.
Output level concerning difference density matrices.
integer > 0 will dump a lot of information—be careful!
Direct SCF procedures build the Fock matrix by exploiting information from previous
iterations for better efficiency. With this keyword information from the last integer
iterations will be used. This feature is switched on with the default value 20, even if the
keyword is absent. The user may reduce the number of iterations employed to smaller values
(e. g. 10) in cases were numerical stability could become an issue. With the value 0
this feature is switched off; the Fock matrix is constructed from scratch in each
iteration.
Control block for convergence acceleration via Pulay’s DIIS
*.
Options are:
specifies the kind of error vector to be used (two different kind of DIIS algorithms)
uses (FDS - SDF) as error vector.
no DIIS
use S-1∕2FDS1∕2 - transposed
Further suboptions:
maximum number of iterations used for extrapolation.
debug level (default: 0)
print applied DIIS coefficients
print DIIS matrix and eigenvalues, too
scaling factor in DIIS procedure: qscal > 1 implies some damping, qscal =
1.0: straight DIIS.
directs the reduction of qscal to qscal = 1.0 (no damping in DIIS), done if
||errvec|| < thrd.
Defaults for $prediag (see above) and $scfdiis
errvec=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.
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).
Disc space specification for two-electron integrals. The following suboptions are available
(and necessary):
Fortran unit number for this file. Unit numbers 30,31,… are recommended.
Filespace in megabytes for this file. size=0 leads to a fully direct run. size is
set by a statistics run, see $statistics. DSCF switches to direct mode if the
file space is exhausted.
Filename. This may also be a complete path name, if you want to store the
integrals in a special directory. Make sure the file is local, otherwise integrals
are transmitted over the network.
Thus your data group $scfintunit may look like this:
Maximum number of SCF iterations (default: 30).
Input/output data group for SCF MOs. You can specify
To perform a calculation without a start vector (i.e. use a core Hamiltonian
guess).
The file where the MOs are written on output (default: mos).
These two options can also be used for $uhfmo_alpha and $uhfmo_beta to use a core guess and write the molecular orbitals to file.
After running define or a TURBOMOLE calculation additional options may appear specifying the origin of the MOs:
These MOs were obtained by projection form another basis set. They should
not be used for wavefunction analysis.
The MOs are converged SCF MOs , the convergence criterion applied was
10-integer
The MOs are unconverged SCF MOs which were written on this data group
after iteration integer. The latter three options are mutually exclusive.
This specifies the FORTRAN format specification which was used for MO
output. The standard format is (4d20.14). (See data group $mo output
format.)
Example:
Your data group $scfmo could look like this after a successful TURBOMOLE run
:
To assist convergence, either the energies of unoccupied MOs can be shifted to higher
energies or, in open-shell cases, the energies of closed-shell MOs to lower energies. In general
a large shift may help to get better convergence.
Options are:
Automatic virtual shell shift switched off.
Automatic virtual shell shift switched on; the energies of virtual orbitals will
be shifted if the HOMO-LUMO gap drops below real such that a gap of real is
sustained. This is the default setting if the keyword is missing with real=0.1.
Option for open-shell cases. Closed shells are shifted to lower energies by real.
The default shift value is closedshell=0.4.
Note: Normally this will disable the automatic shift of energies of virtual
orbitals! To override this, you should append an exclamation mark to the
’automatic’ switch, i.e. specify ’automatic! real’.
Set shifts for special occupied MOs. To use this option, start the line with the
symmetry label and the list of MOs within this symmetry and append the desired
shift in brackets as in the following example:
Integral evaluation threshold. Integrals smaller than real will not be evaluated. Note that
this threshold may affect accuracy and the convergence properties if it is chosen too large. If
$scftol is absent, a default value will be taken obtained from $scfconv by
real = (#bf = number of basisfunctions).
The scratch files allocated by dscf can be placed anywhere in your file systems instead of
the working directory by referencing their pathnames in this data group. All possible
scratch files are listed in the following example:
The following options are allowed
Do not perform integrals statistics
see PARALLEL PROCESSING
Options dscf parallel, grad, mpgrad will be described in the related chapters.
If $statistics dscf has been given integral prescreening will be performed (which is an
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.
Integral storage parameter, which is related to the time needed to calculate the integral.
The larger integer the less integrals will be stored. The default value is integer = 5. (see also
$thize, $statistics)
Integral storage parameter, that determines, together with $thime, the number of integrals
stored on disc. Only integrals larger than real will be stored. The default value is real =
0.100E-04.
Specification of MO occupation for RHF, e.g.
MO occupation of open shells and number of open shells. ’type=1’ here means that there is
only a single open shell consisting e.g. of two MOs:
This data group is necessary for ROHF calculations with more than one open shell.
Example:
For 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:
For further information on ROHF calculations (e.g. with more than one open shell), see the sample input in Section 21.6 and the tables of Roothaan parameters in Section 6.3.
Note that this keyword toggles the ROHF mode also for more than one open shell. If it is not given, the open-shell electrons are simply ignored.
these two data groups specify the occupation of alpha and beta spin UHF MOs (syntax
as any data group related with orbital occupation information, e.g. $closed
shells)
Example:
directs the program to carry out a UHF run. $uhf overwrites closed-shell occupation
specification.
These two data groups contain the UHF MO vectors for alpha and beta spin respectively
(same syntax as $scfmo).
for DFT calculations one has to specify the functional and the grid (for the quadrature of the exchange correlation part). The settings above are default: both lines can be left out if the B-P86 functional and grid m3 are required. Other useful functionals supported are:
(equivalent to the Gaussian98 keyword B3LYP with VWNIII)
(equivalent to the Gaussian98 keyword SVWN with VWNIII)
Possible grids are 1–5 and m3–m5 where grid 1 is coarse (least accurate) and 5 most dense. We recommend however the use of so-called multiple grids m3–m5: SCF iterations with grid 1–3, final energy and gradient with grid 3–5. Usually m3 is fine: for large or delicate systems, try m4. For a reference calculation with a very fine grid and very tight thresholds use ’reference’ as grid specification instead of ’gridsize xy’.
Note: the functionals b3-lyp_Gaussian and s-vwn_Gaussian are made available only for comparability with Gaussian. The functional VWNIII is much less well founded than VWN5 and the TURBOMOLE team does not recommend the use of VWNIII.
Dscf does not run with the keyword $rij: you must call the RI modules Ridft and Rdgrad for energy and gradient calculations. However, it does run with the keyword $rik, but it will ignore all RI settings and do a conventional non-RI Hartree–Fock or DFT calculation.
Enforces an RI-J calculation if module ridft is used, can be used for Hartree-Fock
as well as for DFT calculations with pure or hybrid functionals.
Obsolete keyword - use $rij instead!
Enforces a RI-JK calculation if module ridft is used, can be used for Hartree-Fock
as well as for DFT calculations with pure or hybrid functionals.
Choose the memory core available (in megabyte) for special arrays in the RI
calculation (the less memory you give the more integrals are treated directly,
i.e. recomputed on the fly in every iteration)
Cross reference for the file specifying the auxiliary basis as referenced in $atoms. We
strongly recommend using auxbasis sets optimized for the respective MO basis sets,
e.g. use SVP (or TZVP) for the basis and the corresponding auxbasis as provided by
define (default: file=auxbasis).
Calculation of atomic charges according to the s partial wave and atomic dipole
moments according to the p partial wave as resulting from the auxbasis representation
of the density
If the keyword $rik is found in the control file, ridft performs a Hartree–Fock–SCF calculation using the RI-approximation for both Coulomb and HF-exchange (efficient for large basis sets). For this purpose needed (apart from $ricore):
Cross reference for the file specifying the JK-auxiliary basis as referenced in $atoms.
This group is created by the rijk menu in define.
Multipole-Accelerated-Resolution-of-Identity-J. This method partitions the Coulomb interactions in the near- and far-field parts. The calculation of the far-field part is performed by application of the multipole expansions and the near-field part is evaluated employing the RI-J approximation. It speeds up calculation of the Coulomb term for large systems. It can only be used with the ridft module and requires setting of the $ridft keyword.
The following options are available:
specifies precision parameter for the multipole expansions. Low-precision MARI-J calculations require 1⋅10-6, which is the default. For higher precision calculations it should be set to 1 ⋅ 10-8–1 ⋅ 10-9.
maximum l-moment of multipole expansions. It should be set to a value equal at least twice the maximum angular momentum of basis functions. Default value is 10 and it should probably never be set higher than 18.
Threshold for moment summation. For highly accurate calculations it should be set to 1 ⋅ 10-24.
number of bins per atom for partitioning of electron densities. Default value is 8 and hardly ever needs to be changed.
minimum separation between bins. Only bins separated more than the sum of their extents plus wsindex are considered as far-field. Default is 0.0 and should be changed only by the experts.
maximum extent for charge distributions of partitioned densities. Extents with values larger then this are set to extmax. Hardly ever needs to be changed.
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:
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.
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.
Otherwise the LHF functional can be selected in define: in this case default options are used.
Options for the LHF potential can be specified as follows ( see also lhfprep -help)
calculation of the KLI exchange potential. By default the LHF exchange
potential is computed (off-diag on).
the Slater potential is calculated numerically everywhere: this is more accurate
but much more expensive. When ECP are used, turn on this option.
leads to accurate results only for first-row elements or if an uncontracted basis
set or a basis set with special additional contractions is used: in other cases
numerical-slater on has to be used (this is default).
for asymptotic treatment there are three options:
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.
Full asymptotic-treatment and use of the numerical Slater in the near
asymptotic-region.
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.
the converged Slater and correction potentials for all grid points are saved in the files
slater.pot and corrct.pot, respectively. Using pot-file load, the
Slater potential is not calculated but read from slater.pot (the correction
potential is instead recalculated). For spin unrestricted calculations the
corresponding files are slaterA.pot, slaterB.pot, corrctA.pot and
correctB.pot.
allows the user to specify which occupied orbital will not be included in the
calculation of correction potential: by default the highest occupied orbital is selected.
This option is useful for those systems where the HOMO of the starting orbitals
(e.g. EHT, HF) is different from the final LHF HOMO. homob is for the beta
spin.
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
0…RF - ΔF : basis-set correction potential
RF - ΔF…RF + ΔF : smooth region
RF + ΔF… + ∞ : asymptotic correction
Defaults: RF = 10, ΔF = 0.5
slater-region RNΔNR′FΔ′F
0…RN - ΔN : basis-set Slater potential
RN - ΔN…RN + ΔN : smoothing region
RN + ΔN…R′F - Δ′F : numerical Slater
R′F - Δ′F…R′F + Δ′F : smoothing region
R′F + Δ′F… + ∞ : asymptotic Slater
Note: R′F - Δ′F ≤ RF - ΔF
Defaults: RN = 7, ΔN = 0.5, R′F = 10, Δ′F = 0.5
Use correct-b-region and slater-b-region for the beta spin.
Self-consistent two-component calculations (e.g. for spin-orbit interactions) can be carried out using the module ridft . The following keywords are valid:
enforces two-component-SCF calculations; this option is combinable with $rij, $rik
and $dft.
switches on Kramers-restricted formalism
switches on collinear two-component formalism (not rotational invariant)
enforces DIIS for complex Fock operator.
Relativistic all-electron calculations can be done employing the X2C, the BSS or the DKH Hamiltonian. Implemented for modules dscf and ridft.
switches on X2C calculation.
switches on BSS calculation.
switches on DKH calculation of order integer.
selects parameterization of the DKH Hamiltonian. Valid values are 1 (=default), 2, 3, 4, and
5.
Note in particular that the parametrization does not affect the Hamiltonian up to fourth order. Therefore, as long as you run calculations with DKH Hamiltonians below 5th order you may use any symbol for the parametrization as they would all yield the same results. Higher-order DKH Hamiltonians depend slightly on the chosen paramterization of the unitary transformations applied in order to decouple the Dirac Hamiltonian, but this effect can be neglected. For details on the different parametrizations of the unitary transformations see [235].
switches on local DLU approach.
selects parameterization of the local approximation. Valid values are 0 or 1. For details on
the different parametrizations see [92].
All of these keywords are combinable with $soghf.
Specification of position and magnitude of point charges to be included in the Hamiltonian.
Each point charge is defined in the format
with <x>, <y>, <z> being the coordinates and <q> its charge,e.g.
In addition the following optional arguments may be given:
distance threshold for discarding redundant point charges, default value 10-6.
if given, the selfenergy of the point charge array will will be included in the energy
and the gradient
switches off the check for redundant point charges and the default symmetrization.
This option can significantly speed up the point charge treatment if many of them
are involved - use only if the point charges are well distributed and symmetry is C1,
e.g. when they stem from proper MM runs
print all point charges in the output (default is to print the point charges only if less
than 100 charges given)
The Periodic Electrostatic Embedded Cluster Method (PEECM) functionality provides electronic embedding of a finite, quantum mechanical cluster in a periodic, infinite array of point charges. It is implemented within HF and DFT energy and gradient TURBOMOLE modules: dscf, grad, ridft, rdgrad, and escf. Unlike embedding within a finite set of point charges the PEEC method always yields the correct electrostatic (Madelung) potential independent of the electrostatic moments of the point charges field. It is also significantly faster than the traditional finite point charges embedding.
The basic PEECM settings are defined in the $embed block. It can be redirected to an external file using $embed file=<file_name>.
Following keywords are used for the PEECM calculation setup:
Specifies the number of periodic directions. Allowed values for number are 3 for a bulk three-dimensional system, 2 for a two-dimensional surface slab, and 1 for a one-dimensional system. Default value is 3.
Unit cell parameters in a form of six real values |a|, |b|, |c|, α, β, γ, where |a|, |b|, |c| are lengths of the appropriate cell vectors, α is the angle between vectors b and c, β is the angle between vectors a and c, and γ is the angle between vectors a and b. Default are atomic units and degrees. You can specify unit cell parameters in Å and degrees using cell ang.
Content of the unit cell, where label is the label of the point charge Content of the unit cell, where label is the label of the point charge type and x y z are corresponding Cartesian or fractional crystal coordinates. Defaults are Cartesian coordinates and atomic units. You can specify Cartesian coordinates in Å using content ang and fractional coordinates using content frac. Note that Cartesian coordinates assume that the cell vector a is aligned along the x axis, and the vector b on the xy plane.
Atomic coordinates of the piece of the crystal to be replaced by the QM cluster and surrounding isolation shell (ECPs and explicit point charges), where label is the point charge label and x y z are corresponding Cartesian or fractional crystal coordinates. Defaults are Cartesian coordinates and atomic units. You can specify Cartesian coordinates in Å using cluster ang and fractional coordinates using cluster frac.
Values of point charges (for each atom-type) , where label is the point charge label and charge specifies charge in atomic units.
Values of point charges (for each atom), where label is the point charge label and charge specifies charge in atomic units.
Note that charges and ch_list are mutually exclusive. An integer number n can also be appended to charges or ch_list to set the tolerance for charge neutrality violation to 10-n (default n = 5).
Additionally, the following keywords control the accuracy of PEECM calculation:
Maximum order of the multipole expansions in periodic fast multipole method
(PFMM). Default value is 25.
Electrostatic potential at the lattice points resulting from periodic point charges field
will be output if this keyword is present. Default is not to output.
Well-separateness criterion for PFMM. Default is 3.0.
Minimum accuracy for lattice sums in PFMM. Default is 1.0d-8.
The Conductor-like Screening Model (COSMO) is a continuum solvation model, where the solute
molecule forms a cavity within the dielectric continuum of permittivity epsilon that represents
the solvent. A brief description of the method is given in chapter 17.2. The model is
currently implemented for SCF energy and gradient calculations (dscf/ridft and
grad/rdgrad), MP2 energy calculations (rimp2 and mpgrad) and MP2 gradients (rimp2), and
response calculations with escf. The ricc2 implementation is described in section 17.2.
For simple HF or DFT single point calculations or optimizations with standard settings, we
recommend to add the $cosmo keyword to the control file and to skip the rest of this
section.
Please note: due to improvements in the A matrix and cavity setup the COSMO energies and
gradients may differ from older versions (5.7 and older). The use_old_amat option can be used to
calculate energies (not gradients) using the old cavity algorithm of TURBOMOLE 5.7.
The basic COSMO settings are defined in the $cosmo and the $cosmo_atoms block.
Example with default values:
defines a finite permittivity used for scaling of the screening charges. If the option
ion is added to the same input line, the scaling factor for ions f(ε) = with
x = 0 will be used. Alternatively the x value can be set by adding ion=x with x
as a real value.
skips the COSMO segment statistics run and allocates memory for the given
number of segments.
skips the outlying charge correction.
All other parameters affect the generation of the surface and the construction of the A matrix:
number of basis grid points per atom
(allowed values: i = 10 × 3k × 4l + 2 = 12,32,42,92...)
number of segments per atom
(allowed values: i = 10 × 3k × 4l + 2 = 12,32,42,92...)
distance threshold for A matrix elements (Ångstrom)
distance to outer solvent sphere for cavity construction (Ångstrom)
factor for outer cavity construction in the outlying charge correction
pave intersection seams with segments
leave untidy seams between atoms
amplitude of the cavity de-symmetrization
phase of the cavity de-symmetrization
refractive index used for the calculation of vertical excitations and num.
frequencies (the default 1.3 will be used if not set explicitly)
in case of disjunct cavities only the largest contiguous cavity will be used and the
smaller one(s) neglected. This makes sense if an unwanted inner cavity has been
constructed e.g. in the case of fullerenes. Default is to use all cavities.
If the $cosmo keyword is given without further specifications the default parameter are used (recommended). For the generation of the cavity, COSMO also requires the definition of atomic radii. User defined values can be provided in Ångstrom units in the data group $cosmo_atoms, e.g. for a water molecule:
If this section is missing in the control file, the default values defined in the radii.cosmo file (located in $TURBODIR/parameter) are used. A user defined value supersedes this defaults. $cosmo and $cosmo_atoms can be set interactively with the COSMO input program cosmoprep after the usual generation of the TURBOMOLE input.
The COSMO energies and total charges are listed in the result section. E.g.:
The dielectric energy of the system is already included in the total energy. OC corr denotes the outlying charge correction. The last energy entry gives the total outlying charge corrected energy in the old definition used in TURBOMOLE 5.7 and older versions. The COSMO result file, which contains the segment information, energies, and settings, can be set using: $cosmo_out file= filename.cosmo
Isodensity Cavity: This option can be used in HF/DFT single point calculations only. The $cosmo_isodens section defines the settings for the density based cavity setup (see also chapter 17.2). If the $cosmo_isodens keyword is given without suboptions, a scaled iosodensity cavity with default settings will be created. Possible options are:
activates the density based cavity setup. The default values of nspa and nsph
are changed to 162 and 92, respectively. This values are superseded by the user
defined nspa value of the $cosmo section. By default the scaled density method
is used. The atom type dependent density values are read from the radii.cosmo
file (located in $TURBODIR/parameter).
spacing of the marching tetrahedron grid in Å (default: 0.3Å).
use one isodensity value for all atom types (value in a.u.)
The outlying charge correction will be performed with a radii based outer cavity. Therfore, and for the smoothing of the density changes in the intersection areas of the scaled density method, radii are needed as for the standard COSMO cavity. Please note: The isodensity cavity will be constructed only once at the beginning of the SCF calculation. The density constructed from the initial mos will be used (file mos or alpha/beta in case of unrestricted calculations). Because the mos of an initial guess do not provide a good density for the cavity construction, it is necessary to provide mos of a converged SCF calculation (e.g. a COSMO calculation with standard cavity). We recommend the following three steps: perform a standard COSMO calculation, add the isodensity options afterwards, and start the calculation a second time.
Radii based Isosurface Cavity: The $cosmo_isorad section defines the radii defined isosurface cavity construction. The option uses the algorithm of the isodensity cavity construction but the objective function used depends on the cosmo radii instead of the electron density. The default values of nspa and nsph are changed to 162 and 92, respectively. This values are superseded by the user defined nspa value of the $cosmo section. The resulting surface exhibits smoother intersection seams and the segment areas are less diverse than the ones of the standard radii bases cavity construction. Because gradients are not implemented, the radii based isosurface cavity can be used in single point calculations only.
spacing of the marching tetrahedron grid in Å (default: 0.3Å).
COSMO in MP2 Calculations: The iterative COSMO PTED scheme (see chapter 17.2) can be used with the mp2cosmo script. Options are explained in the help message (mp2cosmo -h). Both MP2 modules rimp2 and mpgrad can be utilized. The control file can be prepared by a normal COSMO SCF input followed by a rimp2 or mpgrad input. The PTE gradients can be switched on by using the
keyword (rimp2 only). Again a normal SCF COSMO input followed by a rimp2 input has to be generated. The $cosmo_correlated keyword forces dscf to keep the COSMO information needed for the following MP2 calculation and toggles on the COSMO gradient contribution.
COSMO in Numerical Frequency Calculations: NumForce can handle two types of COSMO frequency calculations. The first uses the normal relaxed COSMO energy and gradient. It can be performed with a standard dscf or ridft COSMO input without further settings. This is the right method to calculate a Hessian for optimizations. The second type, which uses the approach described in chapter 17.2, is implemented for ridft only. The input is the same as in the first case but Numforce has to be called with the -cosmo option. If no solvent refractive index refind=REAL is given in the $cosmo section of the control file the program uses the default (1.3).
COSMO in vertical excitations and polarizabilities: COSMO is implemented in escf and will be switched on automatically by the COSMO keywords of the underlying SCF calculation. The refractive index, used for the fast term screening of the vertical excitations, needs to be defined in the cosmo section of control file (refind=REAL).
DCCOSMO-RS: The DCOSMO-RS model (see chapter 17.2) has been implemented for restricted and unrestricted DFT and HF energy calculations and gradients (programs: dscf/ridft and grad/rdgrad). In addition to the COSMO settings defined at the beginning of this section, the $dcosmo_rs keyword has to be set.
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 25∘C are implemented in the current TURBOMOLE distribution (see parameter subdirectory):
h2o_25.pot
ethanol_25.pot
methanol_25.pot
thf_25.pot
propanone_25.pot
chcl3_25.pot
ccl4_25.pot
acetonitrile_25.pot
nitromethane_25.pot
dimethylsulfoxide_25.pot
diethylether_25.pot
hexane_25.pot
cyclohexane_25.pot
benzene_25.pot
toluene_25.pot
aniline_25.pot
The DCOSMO-RS energies and total charges are listed in the COSMO section of the output:
The outlying charge correction cannot be defined straight forward like in the normal COSMO model. Therefore, the output shows two corrections that can be added to the Total energy. The first one is the correction on the COSMO level (COSMO) and the second is the difference of the DCOSMO-RS dielectric energy calculated form the corrected and the uncorrected COSMO charges, respectively (DCOSMO-RS). The charges are corrected on the COSMO level only. The Total energy includes the Ediel,RS defined in section 17.2. Additionally the combinatorial contribution at infinute dilution of the COSMO-RS model is given in the output. The use of this energy makes sense if the molecule under consideration is different than the used solvent or not component of the solvent mixture, respectively. To be consistent one should only compare energies containing the same contributions, i.e. same outlying charge correction and with or without combinatorial contribution. Please note: the COSMO-RS contribution of the DCOSMO-RS energy depends on the reference state and the COSMO-RS parameterization (used in the calculation of the chosen COSMO-RS potential). Therefore, the DCOSMO-RS energies should not be used in a comparision with the gas phase energy, i.e. the calculation of solvation energies.
riper shares most of the relevant keywords of the dscf and ridft modules. The $dft data group (see 20.2.7) and auxiliary basis sets defined using the keyword $jbas are always required.
For periodic calculations two additional keywords are necessary:
Specifies the number of periodic directions: n = 3 for a 3D periodic bulk solid, n = 2
for a 2D periodic surface slab and n = 1 for a 1D periodic system. The default value
is 0 for a molecular system.
Specifies the unit cell parameters. The number of cell parameters depends on the
periodicity of the system:
For 3D periodic systems six unit cell parameters |a|, |b|, |c|, α, β and γ need to be provided. Here, |a|, |b| and |c| are lengths of the appropriate cell vectors, α is the angle between vectors b and c, β is the angle between vectors a and c, and γ is the angle between vectors a and b. riper assumes that the cell vectors a and b are aligned along the x axis and on the xy plane, respectively.
For 2D periodic systems three surface cell parameters |a|, |b| and γ have to be provided. Here, |a| and |b| are lengths of the appropriate cell vectors and γ is the angle between a and b. riper assumes that the cell vectors a and b are aligned along the x axis and on the xy plane, respectively.
For 1D periodic systems only one parameter specifying the length of the unit cell has to be provided. riper assumes that periodic direction is along the x axis.
By default, the parameters have to be provided in atomic units and degrees. Alternatively, Å can specified using $cell angs.
Alternatively, lattice vectors can be provided. The number of cell parameters depends
on the periodicity of the system:
For 3D periodic systems three (three-dimensional) lattice vectors need to be provided.
For 2D periodic systems two (two-dimensional) lattice vectors have to be provided. riper assumes that the lattice vectors are aligned on the xy plane.
For 1D periodic systems only one parameter specifying the length of the lattice vector has to be provided. riper assumes that periodic direction is along the x axis.
By default, lattice vectors have to be provided in atomic units. Alternatively, Å can be specified using $lattice angs.
Optionally, for periodic systems a k points mesh can be specified:
nkpoints n1 n2 n3
Specifies components along each reciprocal lattice vector of a Γ centered mesh of k points.
In 3D periodic systems each k point is defined by its components k1, k2 and k3 along the
reciprocal lattice vectors b1, b2 and b3 as
| (20.1) |
For 2D periodic systems k3 = 0. In case of 1D periodicity k3 = 0 and k2 = 0. The three components kj (j = 1,2,3) of k are given as
| (20.2) |
with nj (j = 1,2,3) as integer numbers. For 3D periodic systems n1, n2 and n3 have to be specified. The component n3 can be omitted for 2D periodic systems. In case of 1D periodicity only the n1 has to be specified.
In addition, the options kptlines and recipr within the $kpoints group can be used to obtain band structure plots:
Specifies the number of reciprocal space lines for band structure plots.
Specifies lines in the reciprocal space using three real numbers defining the start point of the line and three real numbers defining its end point. Finally, the number of k points along the line is given as an integer number.
In the following example band energies are calculated along four lines, as specified by the keyword kptlines 4. Each line definition starts in a new line with the keyword recipr, followed by three real numbers defining the start point of the line and three real numbers defining its end point. Finally, the number of k points along the line is given as an integer number. Thus, the first line starts at the point (0.500 0.500 0.500), ends at (0.000 0.000 0.000) and contains 40 k points.
The calculated band structure is written to the file bands.xyz. Each line of the file contains five real numbers: the coordinates k1, k2 and k3 of the k point, its length |k| and the corresponding band energy ϵnk.
Calculation of energy first derivatives with respect to the lattice vectors can be requested using the keyword $optcell. The actual values and gradients of lattice vectors are written to the control file in the data group $gradlatt. Alternatively, another file name can be specified using $gradlatt file=.
Keywords within the section $riper can be used to control precision and parameters of algorithms implemented in riper. If the $riper group is absent, the following default values are used:
The following options are available:
Threshold for integrals neglect and for differential overlap when screening basis functions pairs. Probably never needs to be changed.
Flag for energy calculation only, no gradients.
If set to on charge projection of the auxiliary electron density [111] is performed for molecular systems during calculation of the Coulomb term. The charge projection constraints the charge of auxiliary density exactly to the number of electrons in the system. It is required for periodic systems, otherwise the Coulomb energy would be infinitely large. For molecular systems charge projection leads to a slight increase of the RI fitting error. It may be useful in some cases but we have so far not identified any.
Forces orthonormalization of orbital coefficients every northol SCF iteration.
If set to on full diagonalization of the Coulomb metric matrix [111] is performed and used to solve density fitting equations. When diffuse auxiliary basis functions are used the default Cholesky decomposition of the Coulomb metric matrix may fail due to small negative eigenvalues. In this case the slower method based on a full diagonalization of the metric matrix is necessary.
If pqmatdiag is used pqsingtol sets threshold for neglect of small eigenvalues of the Coulomb metric matrix.
Maximum l-moment of multipole expansions used for calculation of the Coulomb term. The default value hardly ever needs to be changed.
Target number of charge distributions per lowest level box of the octree [107]. The default value hardly ever needs to be changed.
Sets the well-separateness criterion [107]. Octree boxes with centers separated more than sum of their lengths times 0.5×wsicl are considered as well-separated. The default hardly ever needs to be changed. wsicl makes sense only for values ≥ 2.0. For wsicl ≤ 3.0 increasing lmaxmom may be necessary for reasonable accuracy.
Precision parameter used to determine basis function extents [107].
If set to on, an additional acceleration method employing local multipole expansions is used. For periodic systems this leads to a significant speedup of calculations, especially for small unit cells and/or diffuse basis functions. Default value is off and on for molecular and periodic systems, respectively.
For locmult set to on the order of local multipole expansions is increased by locmomadd. The default value probably never needs to be changed.
If set to on the low-memory iterative density fitting (LMIDF) scheme is used for solving the RI equations [110] using the preconditioned conjugate gradient (PCG) method. It is implemented for molecular systems only. Default value is off.
If lpcg is used, lcfmmpcg specifies whether the CFMM is applied for evaluation of the matrix-vector products needed for the PCG solver. Not employing CFMM speeds up the calculations but significantly increases memory demand since the full Coulomb metric matrix has to be stored. Default value is on.
Maximum l-moment of multipole expansions for calculations of Coulomb interactions within the PCG algorithm. It should be set to the same or larger value than lmaxmom.
Sets the threshold parameter controlling accuracy of the PCG solver (see [110] for details). Default value is 1.0 ⋅ 10-9. For lower-precision calculations it can be set to 1.0 ⋅ 10-8 but values larger than 1.0 ⋅ 10-7 are not allowed as these lead to large errors in Coulomb energies and occasionally to SCF convergence problems.
char = at, ss or sp
Specifies the type of preconditioner used in the PCG algorithm. Three
types of preconditioners are implemented and are defined explicitly in
Sec. 7. The sp preconditioner is a default one performing consistently
the best among the preconditioners considered. The at preconditioner
is less efficient in decreasing the number of CG iterations needed
for convergence. However, it has negligible memory requirements
and more favorable scaling properties, albeit with a large prefactor.
The ss preconditioner represents a middle ground between the sp
and at preconditioners both in terms of the efficiency and memory
requirements.
Width of the Gaussian smearing in hartree. See ref. [112] for more information. Note, that [112] uses eV as the unit for the width of the Gaussian smearing.
Specifying desnue along with sigma forces occupancy leading to the number of unpaired electrons equal to desnue.
Keywords within the section $pointvalper can be used to evaluate quantities (electron density/molecular orbitals) for visualization on grid points:
The following options are available:
Specifies output format. Currently .plt, .cub, .xyz and .upt extensions are supported (see 7.3.8).
If present, total (and spin for UHF) density is calculated.
Specifies the number of plotted orbitals. For further details see subsection 7.3.8.
Number of unit cell images n1, n2 and n3 in the periodic directions a, b and c, respectively, for which plot data is generated.
Number of grid points n1, n2 and n3 along each periodic direction. If not specified, value 100 is used for each n.
Specifies the distance real in bohr around the system for which plot grid is generated in aperiodic directions. Default value is 5 bohr.
Number of grid points stored in one octree box during density calculations. Default value is 50. For very large systems or high resolution it may be necessary to increase this parameter to avoid memory allocation problems.
Only valid for plt output format. This format uses orthogonal grids. Therefore, for non-orthogonal unit cells grid data is generated for a rectangular box that contains the supercell (unit cell and its periodic images). By default, the values at grid points outside of the supercell are set to zero. For strongly non-orthogonal systems this may lead to large files. The option full switches off the zeroing of values on grid points outside the supercell.
To calculate a simulated density of states (DOS) keyword $dosper is necessary:
The following options are available:
The width of each Gaussian, default value is 0.01 a.u.
Lower/upper bounds for energy in DOS calculation.
Scaling factor for DOS (total and s-, p-, ... contributions).
Resolution (number of points).
Many of the dscf and ridft keywords are also used by grad and rdgrad.
This keyword and corresponding options are required in gradient calculations only in special circumstances. Just $drvopt is fine, no options needed to compute derivatives of the energy with respect to nuclear coordinates within the method specified: SCF, DFT, RIDFT.
If running a DFT gradient calculation, it is possible to include the derivatives of the quadrature weights, to get more accurate results. In normal cases however those effects are marginal. An exception is numerical calculation of frequencies by Numforce, where it is strongly recommended to use the weight derivatives option. The biggest deviations from the uncorrected results are to be expected if doing gradient calculations for elements heavier than Kr using all electron basis sets and very small grids. To use the weight derivatives option, add
weight derivatives
The option
point charges
in $drvopt switches on the evaluation of derivatives with respect to coordinates of point charges. The gradients are written to the file $point_charge_gradients old gradients will be overwritten.
This module calculates analytically harmonic vibrational frequencies within the HF- or (RI)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.
is the keyword for non-default options of gradient and second derivative calculations.
Possibilities in case of the module aoforce are:
to read a complete Hessian from the input file $hessian and perform only the
frequency analysis
to perform an analysis of normal modes in terms of internal coordinates. Details
about this option and the effect of the printlevel (default is 0) are given in
Section 14. The effect of the keyword only is the same as described above.
fixes the amount of core memory to be used for dynamically allocated arrays (here 50 MiB),
about 70% of available memory should be fine, because $maxcor specifies only the
memory used to store derivatives of density and Fock matrices as well as right
hande side vectors for the CPHF equations. For further details see subsection
20.2.2.
sets the convergence criterion for the 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.
fixes the maximum number of Davidson iterations for the solution of the CPHF-equations to
a value of ten. Normal calculations should not need more than eight iterations, but as a
precaution the default value is 25.
forces the program in case of molecules with C1 symmetry not to use 3N - 6(5) symmetry
adapted but all 3N cartesian nuclear displacement vectors. This option may lead to a
moderate speed-up for molecules notedly larger than 1000 basis functions and 100
atoms.
forces the program not to project out translations and rotations when forming a basis of
symmetry adapted molecular displacements. This option may be needed if a Hessian is
required, that contains translation- and rotation-contributions, e.g. for coupling the system
with low cost methods. Output of the unprojected hessian is done on $nprhessian; format
is the same as for conventional $hessian. Output of the corresponding eigenvalues and
eigenvectors is done analogously on $nprvibrational spectrum and $nprvibrational
normal modes.
causes the program to diagonalize a not mass weighted hessian. Output is on $nprhessian,
$nprvibrational spectrum and $nprvibrational normal modes, because projection of
rotations is not possible in this case.
This keyword allows to trace back the effects of isotopic substitution on vibrational
frequencies. The atom(s) for which isotopic substitution is to be investigated are
specified in subsequent lines of the form (atom index) (mass in special isotope),
e.g.
Sets the number of points for interpolation between the two isotopes compared by the
$isosub option to six. Default value is 21.
Keywords for the treatment of only selected nuclear displacement vectors:
CPHF-iteration is done only for distortions, that are IR active.
CPHF-iteration is done only for distortions, that are Raman active.
Calculation of VCD rotational strenghts after preceding mpshiftrun.
This causes a lowest Hessian eigenvalue search to be performed instead of a complete
force constant calculation. The lowest eigenvalue search consists of the calculation
of a guess-Hessian and macro-iterations to find the solution vector(s) for the lowest
eigenvalue(s). In each macro-iteration the CPHF-equations are solved for the present
search vector(s). $les all 1 means that one lowest eigenvalue for each irrep will be
determined, other numbers of lowest eigenvalues per irrep are admissible too.
Different numbers of lowest eigenvalues for different irreps are requested by e.g.
$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.
Triggers the generation of input files for hotFCHT (program to calculate
Franck-Condon-factors by R. Berger and co-workers). See 14.5.
Save the derivative of the density matrix for subsequent use in the module evib. See
15
Force constant calculations on the DFT level prove to be numerically reliable only with large
integration grids or if one includes the effects of quadrature weights. This is done by default—to
prevent this, insert
no weight derivatives
in $dft.
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
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))
or
or
2. TDA (for HF: CI singles) singlet or triplet or spin-unrestricted or spin-flip excitation energies (HF+RI(DFT))
or
or
or
3. Two-component TDDFT/TDA excitation energies of Kramers-restricted closed-shell systems
or
4. Eigenvalues of singlet or triplet or non-real stability matrices (HF+RI(DFT), RHF)
or
or
5. Static polarizability and rotatory dispersion tensors (HF+(RI)DFT, RHF+UHF)
6. Dynamic polarizability and rotatory dispersion tensors (HF+(RI)DFT, RHF+ UHF)
$scfinstab dynpol unit
list of frequencies
where unit can be eV, nm, rcm; default is a.u. (Hartree). For example, to calculate dynamic polarizabilities at 590 nm and 400 i nm (i is the imaginary unit):
The number and symmetry labels of the excited states to be calculated is controlled by the data group $soes. Example:
will yield the 17 lowest excitations in IRREP b1g, the 23 lowest excitations in IRREP eu, and all excitations in IRREP t2g. Specify $soes alln; to calculate the n first excitations in all IRREPS. If n is not specified, all excitations in all IRREPS will be obtained.
During an escf run, a 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.
The maximum amount of core memory to be allocated for the storage of trial vectors
is restricted to n MB. If the memory needed exceeds the threshold given by $rpacor,
a multiple pass algorithm will be used. However, especially for large cases, this will
increase computation time significantly. The default is 200 MB.
The calculated excitation energies and corresponding oscillator strengths are
appended to a file named ’spectrum’. Possible values of unit are eV, nm and cm-1
or rcm. If no unit is specified, excitation energies are given in a.u.
The calculated excitation energies and corresponding rotatory strengths are appended
to a file named ’cdspectrum’. unit can have the same values as in $spectrum.
Flag for generation of UHF start MOs in a triplet instability calculation. The option
will become effective only if there are triplet instabilities in the totally symmetric
IRREP. The optional real number e specifies the approximate second order energy
change in a.u. (default: 0.1).
Enables calculation of dipole polarizability/rotatory dispersion in the velocity gauge.
Active only for pure DFT (no HF exchange).
Enable calculation of oscillator and rotatory strength sum rules at frequencies
specified by list of frequencies in unit unit (see $scfinstab dynpol). Note that the
sums will be taken only over the states specified in $soes.
the vectors are considered as converged if the Euclidean residual norm is less than
10-n. Larger values of n lead to higher accuracy. The default is a residual norm less
than 10-5.
Sets the upper limit for the number of Davidson Iterations to n. Default is n = 25.
The main keyword that switches on a GW calculation. Provided that the response function is calculated setting this keyword will perform a standard G0W0 calculation with default values for the other flags. There are several options which can be added to the $gw keyword, the syntax is:
With the optional entries:
Default: 1. Number of orbitals to calculate gw for. It is set to the number of occupied orbitals + 5 if set smaller than the number of occupied orbitals.
Default: 0.001. Infinitesimal complex energy shift. Negative value switches to calculating at that value but extrapolating to 0 in linear approximation.
Default: qpenergies.dat. Output filename for the quasiparticle energies.
Default: false (not set). If added as option pure rpa responce function is calculated. If not added, the TDDFT response function is calculated and used to screen the coulomb interaction.
The main keyword that switches on a BSE calculation. Provided that the response function is calculated setting this keyword will perform a standard BSE calculation with default values for the other flags. The following optional entries can be added to the $bse keyword:
Default: Not set. If set, the matrices A and W, respectively, are constructed using KS/HF orbital energies instead of GW quasi-particles energies.
Default: qpenergies.dat. This option defines the file name of GW quasi-particle energies which are read in.
Default: Not set. Compute the auxiliary matrix for the static screened interaction iteratively instead of by Cholesky decomposition.
Default: 1.0d-12. Threshold for convergence in case of the iterative computation of the static screened interaction.
Default: 100. Maximum number of iterations in case of the iterative computation of the static screened interaction.
The keyword $rirpa allows to specify the following options,
Number of frequency integration points, n (default is 60).
HF energy calculation is skipped, (HXX = Hartree + eXact (Fock) eXchange).
Generates profiling output.
Switches on the gradients calculation for RI-RPA.
Computes gradients in the direct RI-MP2 limit.
Manual setting of the number of integral blocks, n, in subroutine rirhs.f. This is for
developers; the default is determined with the $maxcor.
egrad uses the same general keywords as escf and grad, see Sections 20.2.12 and 20.2.15.
The state to be optimized is by default the highest excited state specified in $soes. Note that only one IRREP can be treated at the same time in contrast to escf calculations. When the desired excited state is nearly degenerate with another state of the same symmetry it may be necessary to include higher states in the initial calculation of the excitation energy and vector in order to avoid root flipping. This is accomplished by means of the additional keyword
which explicitly enforces that n-th excited state is optimized. n must not be larger than the number of states specified in $soes.
flag to compute Cartesian non-adiabatic coupling vectors between the excited state of interest and the ground state [236]. This option requires the use of weight derivatives in section dft. It is only implemented for C1 symmetry.
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 |
The data group $maxcor adjusts the maximum size of core memory (n in MB) which
will be allocated during the MP2 run. Recommendation: 3/4 of the actual main
memory at most. For further details see subsection 20.2.2.
Calculation of MP2 gradient is omitted, only MP2 energy is calculated.
The data group $freeze specifies frozen orbitals, in the above syntax by irreducible representations. The symmetry-independent and for standard-applications recommended syntax is
This will freeze the 5 lowest occupied and 2 highest virtual orbitals (alpha and beta count as one in UHF cases). Note that degenerate orbitals count twice (e representations), thrice (t representations) etc. Note: In case of gradient calculations frozen core orbitals are not regarded by mpgrad, moreover freezing of virtual orbitals is generally not supported by mpgrad.
All essential data groups for mpgrad may be generated by the preparation tool mp2prep, apart from $maxcor (see above) these are the following:
specifies the number of loops (or ’passes’) over occupied orbitals, n, performed in the
mpgrad run: the more passes the smaller file space requirements—but CPU time will
go up.
The data group $mointunit specifies:
statistics run (estimation of disc space needed) for the adjustment of the file sizes will
be performed.
calculation of MP2 pair correlation energies.
Note that beside the keywords listed below the outcome of the ricc2 program also depends on the settings of most thresholds that influence the integral screening (e.g. $denconv, $scfconv, $scftol) and for the solution of Z vector equation with 4-index integrals (for relaxed properties and gradients) on the settings for integrals storage in semi-direct SCF runs (i.e. $thime, $thize, $scfintunit). For the explanation of these keywords see Section 20.2.7.
Auxiliary basis set for RI approximation.
Freeze orbitals in the calculation of correlation and excitation energies. For details see
Section 20.2.18.
Print level. The default value is 1.
Specify a directory for large intermediate files (typically three-index coulomb
integrals and similar intermediates), which is different from the directory where the
ricc2 program is started.
The data group $maxcor adjusts the maximum size of core memory which will be
allocated during the ricc2 calculation. $maxcor has a large influence on computation
times. It is recommended to set $maxcor to ca. 75–85% of the available physical core
memory. For further details see subsection 20.2.2.
The calculated excitation energies and corresponding oscillator strengths are
appended to a file named ’spectrum’. Possible values of unit are eV, nm and cm-1
or rcm. If no unit is specified, excitation energies are given in a.u.
The calculated excitation energies and corresponding rotatory strengths are appended
to a file named ’cdspectrum’. unit can have the same values as in $spectrum.
The purpose of this data group is twofold: It activates the Laplace-transformed implementation of SOS-MP2 in the ricc2 module (if the sos option has been specified in $ricc2) and it provides the options to specify the technical details for the numerical Laplace-transformation.
Threshold for the numerical integration used for the Laplace transformation
of orbital energy denominators. The grid points for the numerical integration
are determined such that is the remaining root mean squared error (RMSE) of
the Laplace transformation is < 10-conv. By default the threshold is set to the
value of conv given in $ricc2 (see below).
specifies the ab initio models (methods) for ground and excited states and the most important parameters and thresholds for the solution of the cluster equations, linear response equations or eigenvalue problems. If more than one model is given, the corresponding calculations are performed successively. Note: The CCS ground state energy is identical with the SCF reference energy, CCS excitation energies are identical to CIS excitation energies. The MP2 results is equivalent to the result from the rimp2 module. cis(dinf) denotes the iterative CIS(D) variant CIS(D∞).
Request the calculation of the 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.
If the energy only flag is given after the cis(d) keyword, it is assumed that
only excitation energies are requested. This switches on some shortcuts to avoid
the computation of intermediates needed e.g. for the generation of improved
start vectors for CC2.
If the restart flag is set, the program will try to restart the CC2 calculations
from previous solution vectors on file. If the norestart flag is set no restart will
be done. Default is to do a restart for CC2 if and only if the file CCR0--1--1---0
exists. Note: There is no restart possibility for CCS/CIS or MP2/CIS(D).
If the hard_restart flag is set, the program will try to reuse integrals and
intermediates from a previous calculation. This requires that the restart.cc
file has been kept, which contains check sums and some other informations
needed. The hard_restart flag is switched on by default, if the restart.cc
file is present.
The conv parameter gives the convergence threshold for the ground state energy for the iterative coupled-cluster methods as 10-conv. The default value is taken from the data group $deneps.
The oconv parameter gives an additional threshold for the residual of the
cluster equations (vector function). If this parameter is given, the iterations
for the cluster equations are not stopped before the norm of the residual is
< 10-oconv. By default the threshold is set to oconv =conv-1, or 10× deneps
if no input for conv is given.
If the norm of a vector is smaller than 10-lindep, the vector is assumed to be
zero. This threshold is also used to test if a set of vectors is linear dependent.
The default threshold is 10-15.
gives the maximum number of iterations for the solution of the cluster
equations, eigenvalue problems or response equations (default: 25).
is the maximum number of vectors used in the DIIS procedures for ground
state or excitation energies (default: 10).
the maximum dimension of the reduced space in the solution of linear equations
(default: 100).
print level, by default set to 1 or (if given) the the value of the $printlevel
data group.
Fortran print format used to print several results (in particular one-electron
properties and transition moments) to standard output.
specify wavefunction and electronic state for which a geometry optimization is
intended. For this model the gradient will be calculated and the energy and
gradient will be written onto the data groups $energy and $grad. Required
for geometry optimizations using the jobex script. Note, that in the present
version gradients are only available for ground states at the MP2 and CC2
and for excited states at the CC2 level and not for ROHF based open-shell
calculations. Not set by default. The default model is CC2, the default
electronic state the ground state. To obtain gradients for the lowest excited
state (of those included in the excitation energy calculation, but else of arbitrary
multiplicity and symmetry) the short cut s1 can be used. x is treated as
synonym for the ground state.
the opposite–spin scaling factor cos and the same–spin scaling factor css can
be chosen. If scs is set without further input, the SCS parameters cos=6/5 and
css=1/3 are applied. This keyword can presently only be used in connection
with MP2.
the SOS parameters cos=1.3 and css=0.0 are applied. This keyword can
presently only be used in connection with MP2.
request the calculation of the D1 diagnostic for the ground state wavefunction.
Only needed for MP2 (see above for the alternative input option mp2 d1diag).
For all other correlated methods the D1 diagnostic is evaluated by default
(without significant extra costs).
calculates the second-order corrections to the CCSD(T) energy from the
interference-corrected MP2-F12 (INT-MP2-F12) if $rir12 is switched on. It
can be combined either with the mp2 or the ccsd(t) methods. In the latter case,
the CCSD(T)-INT-F12 energy is printed. The intcorr all keyword writes on
the output all pair energies.
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.)
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.)
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.)
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.)
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.)
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.
char=LCG or R12
The corrfac flag determines which correlation factor is used for the geminal
basis. LCG requires the data group $lcg, which contains the information
regarding exponents and coefficients of the linear combination of Gaussians.
char=off or on
The cabsingles flag determines whether or not the single excitations into the
CABS basis are computed.
The CABS singles are computed in any case if the CABS Fock matrix elements
are computed anyway for the F12 calculation (i.e., for ansatz 2 or r12model B
or comaprox F+K).
char=hf, rohf, boys or pipek
The r12orb flag controls which orbitals are used for the F12 geminal
basis functions. With hf the (semi)-canonical Hartree–Fock orbitals are used
(default). For ROHF-based UMP2 calculations rohf orbitals can be used, which
also implies that the $freeze data group options refer to ROHF rather than
semi-canonical orbitals. For closed-shell species, localised orbitals can be used
with either the Boys or Pipek-Mezey method. For the non-(semi)-canonical
options, the r12orb noinv F12 energy correction is evaluated using active
occupied orbitals transformed to the same basis as the orbitals in the geminal
function.
If no_f12metric is selected the coulomb metric is used in the density fitting
scheme to calculate the four index integrals over the operators f12, f12g12, f122
and f12r12. If f12metric is selected the operator’s own metric is used. The
default for the ricc2 program is no_f12metric, while the pnoccsd program
can only be used with f12metric, where it is therefore the default.
In this data group you have to give additional input for calculations on excited states:
the irreducible representation.
spin multiplicity (1 for singlet, 3 for triplet); default: singlet, not needed for
UHF.
the number of excited states to be calculated within this irrep and for this multiplicity.
the number of roots used in preoptimization steps (default: npre = nexc).
the number of start vectors generated or read from file (default: nstart =
npre).
This flag switches on the calculation of oscillator strengths for excited
state—ground state transitions. Setting the parameter states=all is
mandatory for the calculation of transition properties in the present version.
The operators flag can be followed by a list of operators (see below) for which
the transition properties will be calculated. Default is to compute the oscillator
strengths for all components of the dipole operator.
This flag switches on the calculation of oscillator strengths for excited
state—excited state transitions. Specifying the initial and final states via
istates=all and fstates=all is mandatory for the calculation of transition
properties in the present version. The operators flag can be followed by a list
of operators (see below) for which the transition properties will be calculated.
Default is to compute the oscillator strengths for all components of the dipole
operator.
request the calculation of two-photon transition moments between ground and
excited states. It recognizes the optional suboptions states for specifying
the states for which the two-photon transition moments should be computed,
operators for specifying a list of pairs of one-electron transition operators,
and freq for a list of frequencies for one of the photons. If these suboptions
are specified, two-photon transition moments are computed per default for all
requested singlet states, dipole operators (in the length representation) and
photon energies which are equal 1/2 of the transition energies.
request the calculation of derivatives of transition moments between ground
and excited states. It recognizes the optional suboptions states for specifying
the states for which the two-photon transition moments should be computed,
operators for specifying a list of pairs of one-electron transition operators. For
phosphorescence lifetimes operators should be set to (diplen,soc) and freq
to 0.0d0.
requests the calculation of first-order properties for excited states. For the
states option see spectrum option above; for details for the operators input
see below.
request calculation of the gradient for the total energy of an excited state. If
no state is specified, the gradient will be calculated for the lowest excited state
included in the calculation of excitation energies. The simultaneous calculation
of gradients for several state is possible.
convergence threshold for norm of residual vectors in eigenvalue problems is set to 10-conv. If not given, a default value is used, which is chosen as max(10-conv,10-oconv,10-6), where conv refers to the values given in the data group $ricc2.
convergence threshold used for preoptimization of CC2 eigenvectors is set to
10-preopt (default: 3).
threshold (10-thrdiis) for residual norm below which DIIS extrapolation is
switched on in the modified Davidson algorithm for the non-linear CC2
eigenvalue problem (default: 2).
If the flag leftopt is set the left eigenvectors are computed (default is to
compute the right eigenvectors, for test purposes only).
The bothsides flag enforces the calculation of both, the left and the right
eigenvectors (for test purposes only).
The oldnorm flag switches the program to the old normalization of the
eigenvectors and %T1 and %T2 diagnostics which were identical with those
used in the CC response code of the Dalton program.
In this data group you have to give additional input for the calculation of ground state properties and the solution of response equations:
This flag switches on the calculation of ground state first-order properties (expectation values). The operators flag can be followed by a list of operators (see below) for which the first-order properties will be calculated. Default is to compute the components of the dipole and the quadrupole moment. The option unrelaxed_only suppress the calculation of orbital-relaxed first-order properties, which require solution the CPHF-like Z-vector equations. Default is the calculation of unrelaxed and orbital-relaxed first-order properties. The unrelaxed_only option will be ignored, if the calculation of gradients is requested (see gradient option below and geoopt in data group $ricc2).
requests the calculation of ground state second-order properties as e.g. dipole polarizabilities. The operators flag has to be followed by a comma seperated pair of operators. (If more pairs are needed they have to be given with additional sop commands.) Default is to compute all symmetry-allowed elements of the dipole-dipole polarizability. With the freq flag on can specify a frequency (default is to compute static polarizabilities). The relaxed flag switched from the unrelaxed approach, which is used by default, to the orbital-relaxed approach. Note that the orbital-relaxed approach can not only be used in the static limit (freq=0.0d0). For further restrictions for the computation of second-order properties check Chapter 10.5.
require calculation of geometric gradients. In difference to the geoopt keyword
in the data group $ricc2 this can be used to compute gradients for several
methods within a loop over models; but gradients and energies will not
be written to the data groups $grad and $energy as needed for geometry
optimizations. Note, that in the present version gradients are only available for
MP2 and CC2 and only for a closed-shell RHF reference.
convergence threshold for norm of residual vectors in linear response equations
is set to 10-conv. If not given in the $response data group, a default value is
used, which is chosen as max(10-conv,
10-oconv,10-6), where conv and oconv refer to the values given in the data
group $ricc2.
convergence threshold for the norm of the residual vector in the solution of the
Z vector equations will be set to 10-zconv.
use semi-canonical formulation for the calculation of (transition) one-electron
densities. Switched on by default. The semi-canonical formulation is usually
computationally more efficient than the non-canonical formulation. Exceptions
are systems with many nearly degenerate pairs of occupied orbitals, which have
to be treated in a non-canonical way anyway. (See also explanation for thrsemi
below).
use non-canonical formulation for the calculation of (transition) one-electron
densities. Default is to use the semi-canonical formulation.
the threshold for the selection of nearly degenerate pairs of occupied orbitals
which (if contributing to the density) have to be treated in a non-canonical
fashion will be set to 10-thrsemi. If set to tight the semi-canonical algorithm
will become inefficient, if the threshold is to large the algorithm will become
numerical unstable
threshold for preoptimizating the so-called Z vector (i.e. the lagrangian
multipliers for orbital coefficients) with a preceding RI-CPHF calculation with
the cbas auxiliary basis. The RI-CPHF equations will be converged to a residual
error < 10-zpreopt. Default is zpreopt=4. This preoptimization can reduce
significantly the computational costs for the solution of the Z vector equations
for large basis sets, in particular if they contain diffuse basis functions. For
calculations on large molecules with small or medium sized basis sets the
preoptimization becomes inefficient compared to the large effects of integral
screening for the conventional CPHF equations and should be disabled. This
option is automatically disabled for ricc2 calculations based on foregoing
RI-JK Hartree-Fock calculation.
disable the preoptimization of the Z vector by a preceding RI-CPHF calculation
with the cbas basis set. (Note that the preoptimization is automatically
deactivated if the ricc2 calculation is based on a foregoing RI-JK Hartree-Fock
calculation.)
Common options for keywords in the data groups $ricc2, $response, and $excitations:
input of operator labels for first-order properties, transition moments, etc. Currently
implemented operators/labels are
overlap (charge) operator: the integrals evaluated in the AO basis are ⟨μ|ν⟩
dipole operator in length gauge: ⟨μ|riO|ν⟩ with i = x, y, z; the index
O indicates dependency on the origin (for expectation values of charged
molecules), which in the present version is fixed to (0,0,0)
(all three components, individual components can be specified with the
labels xdiplen, ydiplen, zdiplen).
dipole operator in velocity gauge: ⟨μ|∇i|ν⟩
(all three components, individual components can be specified with the
labels xdipvel, ydipvel, zdipvel).
quadrupole operator ⟨μ|riOrjO|ν⟩
(all six components, individual components can be specified with
the labels xxqudlen, xyqudlen, xzqudlen, yyqudlen, yzqudlen,
zzqudlen).
If all six components are present, the program will automatically give the
electronic second moment tensor (which involves only the electronic
contributions) Mij, the isotropic second moment α = trM and the
anisotropy
|
Furthermore the traceless quadrupole moment
|
(including nuclear contributions) is given.
angular momentum ⟨μ|LiO|ν⟩
(all three components, individual components can be specified with the labels
xangmom, yangmom, zangmom).
electronic force on nuclei ⟨μ||ν⟩, where ZI is the charge of the nucleus I and rI is the position vector of the electron relative to the nucleus (all three components for all nuclei: the labels are xnef001, ynef001, znef001, xnef002, etc. where the number depends on the order in the coord file).
specification of states for which transition moments or first-order properties are to be
calculated. The default is all, i.e. the calculations will be done for all excited states for
which excitation energies have been calculated. Alternatively, one can select a subset of
these listed in parentheses, e.g. states=(ag{3} 1,3-5; b1u{1} 1-3; b2u4). This will
select the triplet ag states no. 1, 3, 4, 5 and the singlet b1u states no. 1, 2, 3 and the singlet
(which is default if no {} is found) b2u state no. 4.
The specification of initial and final states for transition properties between excited states is
mandatory. The syntax is analog to the states option, i.e. either all or a list of of states is
required.
Calculate the double-substitution-based diagnostics D2.
Write MP2/CC2 natural occupation numbers and natural orbitals to a file.
Calculate the error functional δRI for the RI approximation of (ai|bj) integrals
and its gradients with respect to exponents and coefficients of the auxiliary basis set as specified in the data group $cbas. The results are written to $egrad scaled by the factor given with the keyword $cgrad and can be used to optimize auxiliary basis sets for RI-MP2 and RI-CC2 calculations (see Section 1.5).
The ccsdf12 program uses a subset of the data groups of the ricc2 program. For F12 calculations it uses in addition the data groups $rir12, $lcg, $cabs, and $jkbas. See the respective sections for further details.
The ccsdf12 program recognizes in the $ricc2 data group the following options to specify wavefunction methods:
The options mp3, mp4, and cccsd request, respectively, MP3, MP4 and CCSD calculations. The option ccsd(t) request a CCSD calculation with the perturbative triples correction, CCSD(T), and as a side result also the CCSD[T] energy will be printed.
In the data group $rir12 it recognizes ccsdapprox as additional option:
defines the approximation to 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.
Note that beside the keywords listed below the outcome of the pnoccsd program also depends on the settings of most thresholds that influence the integral screening (e.g. $denconv, $scfconv, $scftol). For the explanation of these keywords see Section 20.2.7.
Auxiliary basis set for RI approximation.
Freeze orbitals in the calculation of correlation and excitation energies. For details see
Section 20.2.18.
Print level. The default value is 1.
Specify a directory for large intermediate files (typically three-index coulomb integrals
and similar intermediates), which is different from the directory where the program
is started.
The data group $maxcor adjusts the maximum size of core memory which will be
allocated during the pnoccsd run. $maxcor has a large influence on computation
times. It is recommended to set $maxcor to ca. 75–85% of the available physical core
memory. For further details see subsection 20.2.2.
Only needed for test purposes. It sets the accuracy for the numerical 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.
specifies the ab initio model (method). The current release version is restriced to MP2 which is also the default model.
specifies the localization method; possible choices are boys for Foster-Boys, pm for Pipek-Mezey, ibo for intrinsic bond orbitals (IBOs) and none for canonical (deprecated, only meant for testing) orbitals. Default are Foster-Boys orbitals, IBOs are recommended if seperation of σ- and π-type orbitals is important. (Pipek-Mezey orbitals become very delocal with diffuse basis sets.)
switch between to algorithms for the generation of OSV coefficients. Possible choices are full for an (4) scaling direct diagonalization (cmp. [237]) and davidson for an (3) scaling iterative diagonalization (cmp. [20]). davidson is the default choice and recommended.
the maximal dimension of trial vectors in the iterative OSV generation (prepno davidson). For PNO-MP2 the default is 800. The dimension is bounded by the number of active virtual orbitals and except for small systems a much smaller value as the number of virutals is sufficient. Smaller dimensions increase the performance, but than the iterative scheme might not converge and the program must be restarted with an adjusted dimension. For PNO-MP2-F12 the Default is 4000 since a larger reduced space is required to construct the OSX (cmp. [237]). Usually there is no need to touch this parameter.
specifies the PNO truncation threshold. Default value: 10-7.
specifies the OSV truncation threshold. If no given tolosv is set to 0.1×tolpno, which is the recommended value.
specifies the threshold for selecting orbital and pair-specific auxiliary basis sets for the local RI approximation. If not given tolri is set to 107∕12 ×, which is the recommended value.
specifies the energy threshold for selecting the significant pairs. If not given tolpair is set to (0.1 ×tolpno)2∕3
enables (on) or disables (off) for F12 calculations the use of OPNOs for the occupied orbital spaces in the projectors for the three-electron integrals. Default: on
sets for F12 calculations the truncation thresholds for the complementary auxiliary PNOs for the virtual spaces (CAPNOs) in the projectors for the three-electron integrals. If not specified, default values are calculated from the threshold tolpno.
sets for F12 calculations the truncation thresholds for the orbital specific complementary auxiliary virtuals from which the CAPNOs are generated. If not specified, default values are chosen as 0.1 × the tolcapno thresholds, which is the recommended choice.
sets for F12 calculations the truncation thresholds for the selection of PNOs for the occupied orbital spaces (OPNOs) in the projectors for the three-electron integrals. If not specified, default values are calculated from the threshold tolpno.
set for F12 calculations the truncation thresholds for the orbital specific auxiliary occupied orbitals from which the OPNOs are generated. If not specified, default values are chosen as 0.1 × the tolopno threshold, which is the recommended choice.
sets for F12 calculations the pair specific projector AQij for the three-electron
integrals. Possible choices are:
projector 1 for 1 -12 -12 - O12′′-1′′O2
projector 2 for 1 -12 -12 -1V 2′′- V 1′′2
projector 3 for 1 -12 -12 -12′′-1′′2
With projector 3 AQij becomes independend of the system size and is
therefore the default. But note that for the (4) scaling algorithm this can
not be exploited and here projector 1 is more suited.
The conv parameter gives the convergence threshold for the ground state energy as 10-conv. The default threshold is 10-7.
The oconv parameter gives an additional threshold for the residual of the ground state equations as < 10-oconv. The default threshold is 10-3.
If the norm of a vector is smaller than 10-lindep, the vector is assumed to be zero. This threshold is also used to test if a set of pre-PNOs is linear dependent. The default threshold is 10-12.
gives the maximum number of iterations for the solution of the cluster equations, eigenvalue problems or response equations (default: 25).
is the maximum number of vectors used in the DIIS procedures for the ground state equations (default: 10).
not used in the current release.
the opposite–spin scaling factor cos and the same–spin scaling factor css can be chosen. If scs is set without further input, the SCS parameters cos=6/5 and css=1/3 are applied.
the SOS parameters cos=1.3 and css=0.0 are applied.
for the description of this data group see Sec. 20.2.19.
define what kind of nonlinear parameters are to be optimized by relax and specify
some control variables for parameter update.
Available options are:
optimize molecular structures in the space of internal coordinates using
definitions of internal coordinates given in $intdef as described in Section 4.1
( default: on).
optimize molecular structures in redundant internal coordinates using
definitions of redundant internal coordinates given in $redundant. For an
optimization in redundant internal coordinates option internal has to be
switched on too, and option cartesian has to be switched off (default: on).
optimize molecular structures in the space of (symmetry-distinct) cartesian
coordinates (default: off).
optimize basis set exponents (default=off).
Available suboptions are:
exponents of uncontracted basis functions will be optimized after
conversion into their logarithms (this improves the condition of the
approximate force constant matrix obtained by variable metric methods
and the behavior of the optimization procedure); scale factors of
contracted basis functions will not be affected by the logarithm suboption
ALL basis set exponents will be optimized as scale factors
(i.e. contracted blocks and single functions will be treated in the same
way); if both suboptions (scale and logarithm) are given the logarithms
of the scale factors will be optimized
optimize a global scaling factor for all basis set exponents (default: off).
define some variables controlling the update of coordinates.
Available options are:
maximum allowed total change for update of coordinates. The maximum change
of individual coordinate will be limited to dqmax∕2 and the collective change
dq will be damped by dqmax∕⟨dq∣dq⟩ if ⟨dq∣dq⟩ > dqmaxq
(default: 0.3)
calculate geometry update by inter/extrapolation of geometries of the last two
cycles (the interpolate option is always switched on by default, but it is only
active ANY time if steepest descent update has been chosen, i.e. $forceupdate
method=none; otherwise it will only be activated if the DIIS update for the
geometry is expected to fail)
provide a statistics output in each optimization cycle by displaying all (the last
integer, default setting by define is 5) subsequent coordinates, gradient and
energy values (default: on).
the presence of this keyword forces relax to provide informational output about the usage
of DIIS for the update of the molecular geometry.
special input related to the transformation of atomic coordinates between cartesian and
internal coordinate spaces (default: off).
Available options are:
maximum number of iterations for the iterative conversion procedure internal
→ cartesian coordinates (default: 25).
convergence criterion for the coordinate conversion (default: 1.d-10).
this switch activates special tasks: transform coordinates/gradients/ hessians between
spaces of internal/cartesian coordinates using the definitions of internal coordinates
given in $intdef:
available suboptions are:
the direction of the transformation is indicated by the direction of the arrow
this data group defines both the method for updating the approximate force constant
matrix and some control variables needed for the force constant update.
Options for method:
no update (steepest descent)
Murtagh–Sargent update
Davidon–Fletcher–Powell update
Broyden–Fletcher–Goldfarb–Shanno update
combined (bfgs+dfp) update
Schlegel update
Ahlrichs update (macro option)
number of structures used
maximum number of geometries (= rank of the update procedure, for ahlrichs only)
minimum number of geometries needed to start
update
for an explanation see suboptions pulay given below e.g. ahlrichs numgeo=7 mingeo=3 maxgeo=4 modus=<g|dg> dynamic
NOTES:
try to find an optimal linear combination of the coordinates of the numpul previous
optimization cycles as specified by modus (see below).
Available suboptions are:
char=<g|g> or <g|dq> or <dq|dq> defines the quantity to be minimized (dq
= internal coordinate change).
fmode specifies the force constants to be used (only if char=<g|dq> or
<dq|dq>!)
fmode=static: use static force constants
fmode=dynamic: use updated force constants
real defines the threshold for the quantity g * dq∕|g|*|dq| which defines the
angle between gradient vector and coordinate change (default: 0.1). If pulay is
used in connection with a multidimensional BFGS update for the hessian than
the default is real=0.0. If > -real the pulay update for the geometry
is expected to fail and will be ignored. For example:
pulay numpul=4 maxpul=4 minpul=3 modus=<dq|dq> static fail=0.2
options for $forceupdate
update only the diagonal force constants (update for off-diagonals will be
suppressed) (only active if method=ms, dfp, bfgs)
this allows to damp off-diagonal force constants by 1/real (compare offreset,
which discards off-diagonals completely). Only values > 1.0 will be accepted.
This option is active only within one relax run and will be disabled
automatically by relax. This is useful in difficult cases, where the non-diagonal
update has lead to too large non-diagonal elements of the hessian.
reset off-diagonal force constants to zero. This option will be active for the
current optimization cycle only, i.e. it will be removed by relax after having
discarded off-diagonals!
optimization cycle specification of a maximum energy change allowed (given in
mHartree) which will be accepted using the actual approximate force constant
matrix from $forceapprox; if this energy change will be exceeded, the force
constants will be scaled appropriately
(The default: 0.0 means NO action)
scaling factor for the input hessian (default: 1.0).
lower bound for eigenvalues of the approximate hessian (default: 0.005); if any
eigenvalue drops below threig, it will be shifted to a reasonable value defined
by:
default: texttt0.005.
upper bound for eigenvalues of the hessian; if any eigenvalue exceeds thrbig,
it will limited to this value (default: 1000.0).
damp the variable metric update for the hessian by 1∕(1+ real) (default: 0.0).
specify initialization of the (approximate) force constant matrix.
Available options are:
this activates or deactivates initialization; if on has been set, relax will
provide an initial force constant matrix as specified by one of the possible
initialization options as described below and will store this matrix in data group
$forceapprox; after initialization relax resets $forceinit to off!
provide a diagonal force constant matrix with:
available suboptions are:
this will lead to an assignment of diagonal elements (default: 1.0)).
this will lead to an assignment of initial force constant diagonals
depending on the coordinate type.
Provide individual defined force constant diagonals for
This does not work for basis set optimization. For the correct syntax of ‘fdiag=..’ see descriptions of $intdef, $global
read a cartesian (e.g. analytical) hessian from $hessian and use it as a start
force constant matrix; if $optimize internal has been set: use its transform
in internal coordinate space. If large molecules are to be optimized, it may be
necessary (large core memory requirements!) to deactivate the numerical
evaluation of the derivative of the B-matrix with respect to cartesian
coordinates, which is needed to transform H(cart) →H(int) exactly by
specifying no dbdx.
These keywords depend on the optimization task to be processed and are updated by the
corresponding program (i. g. SCF energy).
This data block contains 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:
atomic index followed by diagonal elements of the m-matrix for this atom
The scratch file ftmp allocated by relax can be placed anywhere in your file systems
instead of the working directory by referencing its pathname in this data group as follows:
Definitions of internal coordinates and, optionally, values of internal coordinates
(val=..., given in a.u. or degrees) or force constant diagonal elements (fdiag=...).
Cartesian coordinates and gradients calculated in subsequent optimization cycles.
Entries are accumulated by one of the gradient programs (grad, mpgrad, rimp2, ricc2,
egrad, etc.).
Basis set exponents scale factors and their gradients as calculated in subsequent
optimization cycles. Entries are accumulated by one of the gradient programs.
Global scale factors and gradients as calculated in subsequent optimization cycles.
Entries are accumulated by the grad or aoforce program.
Allows to augment internal SCF gradients by approximate increments obtained
from treatments (e.g. correlation or relativistic) on higher level. See the example
below.
Approximate force constant matrix (as needed for geometry optimization tasks). The
storage format may be specified by the
available options:
the default format is format=(8f10.5), but other 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:
this data block contains the analytical cartesian force constant matrix (with translational
and rotational combinations projected out) as output by the aoforce program and may be
used to supply a high quality force constant matrix $forceapprox for geometry
optimizations (specifying $forceinit on carthess, or $interconversion cartesian –>
internal hessian).
either updated cartesian coordinates if a successful coordinate update has been
performed, or cartesian coordinates for input internal coordinates if only a conversion
from internal to cartesian coordinates has been performed.
updated basis set exponents, basis sets contraction coefficients or scaling factors, if
$optimize basis on has been specified.
updated global scaling factor for all basis set exponents, if $optimize global on has
been specified.
an approximate force constant matrix to be used in quasi-Newton type geometry
optimizations; this matrix will be improved in subsequent optimization cycles if one of
the variable-metric methods ($forceupdate) has been chosen. See 5.3.13 and 20.2.22.
a static (i.e. never updated) approximate force constant matrix to be used
in DIIS-type geometry optimizations. It will be initialized by relax specifying:
$forceupdate pulay …modus=<dq|dq> static.
The next data groups are output by relax (depending on the optimization subject) in order to control the convergence of optimization procedures driven by the shell script jobex.
real is the absolute value of the maximum component of the corresponding gradient.
In order to save the effort for conversion of accumulated geometry and gradient data (as needed for the force constant update or the DIIS update of the geometry) to the optimization space, within which the geometry has to be optimized, one may specify the keyword
Then the relax program accumulates all subsequent coordinates and gradient as
used in optimization in this data group (or a referenced file). This overrides the input
of old coordinate and gradient data from data blocks $grad, $egrad, …as accumulated
by the grad program.
degrees
Only non-default values are written in the control file except:
Following options are available:
Index of the Hessian eigenvector to follow for transition structure search (transition
vector). Eigenpairs are sorted in ascending order, i.e. with increasing eigenvalues and
start with index 1. The eigenpairs corresponding to translations and rotations are
shifted to the end. For minimization the value 0 has to be specified.
Method of hessian update. For minimization default is BFGS, for TS search default
is Powell and none is for no update.
Recompute the full Hessian every N’th step during a transition state search. The
default is zero and the Hessian is read in or computed in the first step only. If the
standard Hessian update methods fail, it can help to use this keyword. Warning: This
will make the calculation much more time demanding!
Freezing transition vector index.
diagonal hessian elements for diagonal Hessian guess (default: 0.5).
Maximum allowed value for trust radius (default: 0.3).
Minimum allowed value for trust radius (default: 1.0d-4).
Initial value for trust radius (default tradius: radmax = 0.3).
threshold for energy change (default: 1.0d-6).
threshold for maximal displacement element (default: 1.0d-3).
threshold for maximal gradient element (default: 1.0d-3).
threshold for RMS of displacement (RMS = root mean square)
(default = 5.0d-4)
threshold for RMS of gradient (default: 5.0d-4).
All values are in atomic units.
$properties specifies the global tasks for program moloch by virtue of the following options
a missing option or a option followed by the flag off will not be taken into account. The flag active may be omitted. For most of these options (with the only exceptions of trace and cowan-griffin), there are additional data groups allowing for more detailed specifications, as explained below.
if moment is active you need
to compute the 0th, 1st, 2nd and 3rd moment at the reference point 0 0 0.
if potential is active you need
to compute the electrostatic potential (pot) and/or electrostatic field (fld) and/or electrostatic field gradient (fldgrd) and/or the zeroth order contribution to the diamagnetic shielding (shld) at reference point 0 0 0.
if localization is active you need $boys to perform a 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:
if population analyses is active you need
to perform a Mulliken population analysis. The options specify the output data:
print molecular orbital contributions to atomic s, p, d,…-populations
print molecular orbital contributions to overlap populations
print atomic netto populations
print contributions of (irreducible) representations to atomic s,p,d,…-populations
print contributions of (irreducible) representations to overlap populations
or
$loewdin
to perform a Löwdin population analysis (options are invalid here). A Löwdin population
analysis is based on decomposing D instead of DS in case of a Mulliken
PA.
or
to perform a population analysis based on occupation numbers (the options are not necessary and produce some output data concerning the modified atomic orbitals):
print MO contributions to occupation numbers of modified atomic orbitals (MAOs).
print all MAOs on standard output
print all MAOs to file mao.
This kind of population analysis basically aims at 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:
atom symbol
list of all atoms for which this definition should apply. The syntax for this list is as usual in TURBOMOLE, e.g. 2,3,8-10,12
means number of MAOs to be included
means selection criterion for MAOs. meth can be occ (default), eig, or man
string, where occ denotes selection of MAOs by occupation numbers, eig
selection by eigenvalues and man allows manual selection. In the latter case
the string (max. 8 characters) appended to man serves as nickname for the
definition of the MAOs to be chosen. This nickname is expected to appear
as the leftmost word in a line somewhere within data group $mao selection
and is followed by the indices of the modified atomic orbitals which are to be
selected.
means the threshold to be applied for the selection criteria occ or eig (default:
0.1).
Example:
option plot is out of fashion; to plot quantities on a grid, rather use $pointval in
connection with dscf, ridft, rimp2 or egrad, as described below. If nevertheless plot is
active you need
to obtain two-dimensional plot data of mo 4a1g (the plane is specified by origin and two vectors with grid range and number of grid points) which is written to file 4a1g. Several plots may be obtained (#1, #2 etc.) at the same time. Use tool ’konto’ to visualize the plot.
Note: This is the old-fashioned way to plot MOs and densities. A new—and easier—one is to use $pointval, as described below.
if fit is active you need
Each line refers to all atoms, the line specifies a spherical layer of grid points around the atoms. The number of points and their distance from the van der Waals surface [Bohr] are given (the default is 1.0).
one line only, smoothing of the layers of grid points around the molecule: the real number is used to define isopotential surfaces on which the points of the layers have to lie.
One line per element has to be specified, it contains the name of the element and the van der Waals radius in [Bohr].
Properties of RHF, UHF and (two-component) GHF wave functions as well as those of SCF+MP2 densities or such from excited state DFT-calculations can be directly analyzed within the respective programs (dscf, ridft, mpgrad, rimp2 and egrad). In case of spin-unrestricted calculations results are given for total densities (Dα + Dβ) and spin densities (Dα - Dβ). If not explicitly noted otherwise, in the following "D" is the SCF density, D(SCF), in case of dscf and ridft, the MP2-corrected density, D(SCF)+D(MP2), for mpgrad and rimp2 and the entire density of the excited state in case of egrad. For modules dscf and ridft the analysis of properties may be directly started by calling dscf -proper (or ridft -proper). In case of mpgrad and rimp2 this is possible only, if the MP2 density has already been generated, i.e. after a complete run of mpgrad or rimp2.
Functionalities of analyses are driven by the following keywords.
leads to calculation of relativistic corrections for the SCF total density in case of
dscf and ridft, for the SCF+MP2 density in case of rimp2 and mpgrad and for that
of the calculated excited state in case of egrad. Quantities calculated are expectation
values < p2 >,< p4 > and the Darwin term (∑
1∕ZA * ρ(RA)).
yields calculation of electrostatic moments arising from nuclear charges and total
electron densities. Also without setting this keyword moments up to quadrupole are
calculated, with respect to reference point (0,0,0). Supported extensions:
By integer i; the maximum order of moments is specified, maximum and default is i=3 (octopole moments), real numbers x,y,z allow for the specification of one or more reference points.
drives the options for population analyses. By default a Mulliken PA in the basis of
cartesian atomic orbitals (CAOs) is performed for the total density (Dα + Dβ) leading to
Mulliken (brutto) charges and, in case of spin-unrestricted calculations also for
the spin density (Dα - Dβ) leading to Mulliken (brutto) numbers for unpaired
electrons. Besides total numbers also contributions from s-, p-, …functions are listed
separately.
Two-component wavefunctions (only module ridft and only if $soghf is set): In two-component calculations instead of Sz |(Sx,Sy,Sz)| is written to the output. Additionally a vector-file named spinvec.txt is written, which includes the resulting spinvector for each atom in the system (also the direction).
The following modifications and extensions are supported, if the respective commands are written in the same line as $pop:
Additional information about px,py,pz (and analogous for d and f functions) is displayed (lengthy output).
Contributions are plotted only if arising from atoms selected by list.
Contributions smaller than thrpl are not displayed (default: 0.01).
Mulliken atomic overlap matrix is displayed.
Mulliken netto populations (diagonal elements of Mulliken overlap matrix) are calculated.
Summed Mulliken contributions for a group of molecular orbitals defined by
numbers referring to the numbering obtained e.g. from the tool eiger. Note
that occupancy of MOs is ignored, i.e. all orbitals are treated as occupied.
Mulliken contributions for single MOs defined by numbers (independent of whether
they are occupied or not). If this option is valid, one may additionally
set
to calculate a (simulated) density of states by broadening the discrete
energy levels with Gaussians and superimposing them. The width of
each Gaussian may be set by input (default: 0.01 a.u.). The resolution
(number of points) may be chosen automatically (default values are
usually sufficient to generate a satisfactory plot) or specified by hand.
The output files (dos in case of RHF wave functions, and dos_a+b,
dos_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:
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.
to perform a natural population analyses [213]. The possible options (specified in the same
line) are
AO | must be provided, the CAO case is not implemented. |
tw=real | Threshold tw to circumvent numerical difficulties in computing Ow |
(default: tw=1.d-6). |
idbgl=integer | Debug level |
(default: idbgl=0). |
ab | For UHF cases: Print alpha and beta density results. |
short | Print only natural electron configuration and summary. |
Example:
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:
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.
to perform a population analyses based on occupation numbers [214] yielding
"shared electron numbers (SENs)" and multicenter contributions. For this method
always the total density is used, i.e. the sum of alpha and beta densities in case
of UHF, the SCF+MP2-density in case of MP2 and the GHF total density for
(two-component-)GHF.
The results of such an analysis may depend on the choice of the number of modified atomic orbitals ("MAOs"), which can be specified by an additional line; without further specification their number is calculated by the method "mix", see below. Note: One should carefully read the information concerning MAOs given in the output before looking at the numbers for atomic charges and shared electron numbers.
to specify how MAOs are selected per atom.
Available options are:
a) for the way of sorting MAOs of each atom:
MAOs are sorted according to their eigenvalue (those with largest EW
finally are chosen). This is the default.
MAOs are sorted according to their occupation; note that the number
of all occupation is NOT the number of electrons in the system. This
option is kept rather for historical reasons.
b) for the determination of the number of MAOs:
A fixed number of MAOs is taken for each atom; usually this is the
number of shells up to the complete valence shell, e.g. 5 for B-Ne, 6
for Na-Mg, etc. Exceptions are Elements Sc (Y, La), Ti (Zr, Hf), V
(Nb, Ta) for which not all five d-shells are included, but only 2, 3 or
4, respectively. This modification leads to better agreement with partial
charges calculated by an ESP-fit.
All MAOs with an eigenvalue larger than <real> are chosen; default
is <real>=0.1. This and the following two options are not valid in
connection with occ.
Maximum of numbers calculated from fix and thr=0.1 is taken.
2:1 mixture of fix and thr=0.1. This choice gives best agreement
(statistical) with charges from an ESP-fit. It is the default choice.
c) for additional information about MAOs:
Eigenvalues and occupations for each MAO are written to output.
Entire information about each MAO is written to output. Lengthy.
Further for each atom the number of MAOs and the sorting mode can be specified individually in lines below this keyword. Example:
atom 1,3-4 eig 7
leads to choice of the 7 MAOs with largest eigenvalue at atoms 1, 3-4.
enables the generation of localized molecular orbitals (LMOs) using Boys localization. By
default, all occupied orbitals are included, localized orbitals are written (by default
in the AO basis) to file(s) lmo in case of RHF and lalp and lbet in case of
UHF orbitals. Note, that LMOs usually break the molecular symmetry; so, even
for symmetric cases the AO (not the SAO) basis is used for the output. The
localized orbitals are sorted with respect to the corresponding diagonal element
of the Fock matrix in the LMO basis. In order to characterize these orbitals,
dominant contributions of (canonical) MOs are written to standard output as
well as results of a Mulliken PA for each LMO (for plotting of LMOs see option
$pointval).
The keyword allows for following options (to be written in the same line):
Include only selected MOs (e.g. valence MOs) in localization procedure
(numbering as available from Eiger).
maximum number of orbital rotations to get LMOs; default value is 10000
(sometimes not enough, in particular for highly delocalised systems).
lower threshold for displaying MO and Mulliken contributions (default: 0.1).
LMOs are written to file in the CAO basis (instead of AO).
Pipek-Mezey localization is used.
Foster-Boys localization is used.
Intrinsic Bond Orbitals are used.
Compute LMO centers, i.e. the diagonal matrix elements of the position
operator i = ⟨i||i⟩
Compute LMO spreads defined as σi = = ,
where i is the LMO center
sort LMOs such that LMOs with spatially close centers are also close in
index space, by default the LMOs are sorted according to the orbital energy
expectation values
print timings for localization procedures
triggers the generations of a wfn file. It can be used in dscf/ridft single-point calculations or
in ricc2/egrad gradient calculations.
fits point charges at the positions of nuclei to electrostatic potential arising from electric
charge distribution (also possible for two-component calculations, for UHF cases also for
spin density). For this purpose the ("real") electrostatic potential is calculated at spherical
shells of grid points around the atoms. By default, Bragg-Slater radii, rBS, are taken as
shell radii, for each atom the number of points is given by 1000 ⋅ rBS2, the total
number of points is the sum of points for each atom reduced by the number of
points of overlapping spheres. Non-default shells (one or more) can be specified as
follows:
$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.
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
By default the origin is at (0,0,0), the radius is chosen large enough to include the whole 3D box and all contributions are regarded (eps=0).
Data different from total and spin densities are generated by following (combinable) settings (to be written in the same line as statement $pointval):
leads to calculation of electrostatic potential arising from electron densities, nuclei and—if present—constant electric fields and point charges. The densities used for calculation of potentials are the same as above; the respective filenames are generated from those of densities by replacement of the "d" (for density) by a "p" (for potential). By "pot eonly" only the electronic contribution to the electrostatic potential is calculated.
calculation of electric field. Note, that for 3D default output format (.plt, see below) only norm is displayed. Densities used are the same as above, filenames are generated from those of densities by replacement of "d" (for density) by "f" (for field).
calculation of amplitudes of MOs specified by numbers referring to the
numbering obtained e.g. from the tool eiger in the same format. The
respective filenames are 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.
calculation of amplitudes of LMOs (previously generated by $localize)
ordered by the corresponding diagonal element of the Fock matrix in the LMO
basis.
calculation of amplitudes of NMOs (previously generated by
$natural orbitals file=natural and
$natural orbital occupation file=natural
has to be set, if additionally to one of the above quantities also the density is to be computed.
calculation of the Kohn-Sham exchange-correlation potential. It is only valid for DFT calculations and it works for all exchange-correlation functionals, including LHF. Note that for hybrid functionals, only the Kohn-Sham part of the potential will be computed (the HF part is a non-local-operator and can’t be plotted). For GGA functional the full potential will be computed (local and non-local terms)
For line plots the output file is tx.vec . For UHF calculations the output files are tx.vec (alpha-spin potential) and sx.vec (beta-spin potential).
For a line plot the file has three columns: 1: total potential 2: local term ( or Slater-potential for LHF) 3: non-local terms or Correction term for LHF
Output formats may be specified by e.g. fmt=xyz if written to the same line as $pointval. Supported are:
in case of scalars (density, (L)MO amplitudes, electrostatic potential) this format is: (x,y,z,f(x,y,z)). In case of vectors components of the vector and its norm are displayed. This format is valid for all types of grid (3D, plane, line, points, see below), it is the default format in case of calculation of values at single points. Output file suffix is .xyz.
only for 3D, default in this case. Data are written to binary files that can be directly read by gOpenMol. Note, that this output is restricted to scalar quantities; thus in case of vectors (E-field) only the norm is plotted. Output file suffix is .plt.
only for 3D. Data are written to ASCII files that can be imported by e.g. gOpenMol. Note, that this output is restricted to scalar quantities; thus in case of vectors (E-field) only the norm is plotted. Output file suffix is .map.
a format compatible with gOpenMol for visualization of vectors v. The format is x,y,z,vx,vy,vz.
a format which writes out all data in extended precision as ASCII file. Useful to import in general 3D visualization tools or spread sheet programs. Self-explaining format (open one in an editor to see its definition).
for planes and lines (default in these cases). In case of a line specified by α ⋅ (see below) output is α,f(x,y,z) for scalars, for vectors components and norm are displayed. vectors. Analogously, in case of planes it is α,β,f(x,y,z). The output (file suffix .vec) may be visualized by plotting programs suited for two-dimensional plots. A command file (termed gnuset) to get a contour plot by gnuplot is automatically generated.
only for 3D, writes out data in Cube format which can be imported by many external visualization programs.
For 3D grids non-default boundarys, basis vector directions, origin and resolutions may be specified as follows:
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:
This will shift the origin of the plot into the center of mass of the specified set of atoms and align the grid axes with their principal axes: If two or more atoms are specified which lie on one line grid axis 1 is rotated into this line. If the specified subset consists of three atoms located in a plane the grid axis 1 and 2 are rotated into this plane. If more than three atoms are specified the grid axis 1 and 2 are rotated into a plane which fitted to the position of all specified atoms. With the rotate option the grid can be rotated around a grid axis. It accepts as input the rotation angle the index (1/2/3) of the grid axis around which the grid will be rotated.
Grids of lower dimensionality may be specified (in the same line as $pointval) by typing either geo=plane or geo=line or geo=point The way to use is best explained by some examples:
Values are calculated at a plane spanned by vectors (0,1,0) and (0,0,1) centered at (1,1,1).
Values are calculated at a line in direction (0,1,0) centered at (0,0,1). Output format as above.
Values are calculated at the two points (7.0,5.0,3.0) and (0.0,0.0,7.0).
Plane-averaged density can be computed by
The generated file td.vec will contain the quantity
| (20.3) |
The ab initio molecular dynamics (MD) program frog needs a command file named mdmaster. The interactive Mdprep program manages the generation of mdmaster and associated files. It is always a good idea to let Mdprep check over mdmaster before starting an MD run. Mdprep has 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:
Number of MD time steps to be carried out. $nsteps is decreased by 1 every time
frog is run and JOBEX -md stops when $nsteps reaches 0.
Number of atoms in system.
The file containing the current position, velocity, time and timestep, that is, the input
configuration. During an MD run the $current information is generally kept at the
end of the $log file.
The file to which the trajectory should be logged, i.e. the output: t=time (a.u.);
atomic positions x,y,z (Bohr) and symbols at t;
timestep (au) Δt;
atomic symbols and velocities x,y,z (au) at t - (Δt∕2);
kinetic energy (H) interpolated at t, ab initio potential energy (H) calculated at t,
and pressure recorded at the barrier surface (atomic units, 1 au = 29.421 TPa) during
the corresponding timestep;
ab initio potential energy gradients x,y,z (H/Bohr) at t.
This file can be manipulated with log2? tools after the MD run (Section 1.5).
The status of the MD run is a record of the action carried out during the previous MD
step, along with the duration of that step. The format matches that of $md_action
below.
Canonical dynamics is supported using the Nosé-Hoover thermostat. This option can be enabled in Mdprep or by the following syntax:
Here, T specifies the temperature of the thermostat in K (500 K in the example) and t specifies the thermostat relaxation time in a.u. (100 a.u. in the example). It is advisable to choose the thermostat relaxation 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:
Integer random number seed
Arbitrary title
To determine the trends in kinetic energy and total energy (average values and overall drifts) it is necessary to read the history of energy statistics over the recent MD steps. The number of MD steps recorded so far in each log file are therefore kept in the $log_history entry: this is updated by the program each step. The length of records needed for reliable statistics and the number of steps over which changes are made to kinetic energy (response) are specified in $ke_control.
$barrier specifies a virtual cavity for simulating condensed phases. The optional flag, angstroms, can be used to indicate that data will be entered in Ångstrøms rather than Bohr.
can be one of orth, elps, or none, for orthorhombic, ellipsoidal, or no barrier
(the default) respectively.
are the +x,y,z sizes of the cavity. In this case, an ellipsoid with a major axis of
20 Å along y, semi-major of 15 Å on z, and minor of 10 Å on x.
is the Hooke’s Law force constant in atomic units of force (H/Bohr) per length
unit. Here, it is 2.0 H/Bohr/Ångstrøm, a bastard combination of units.
is the effective limit to the restorative force of the barrier. For this system, an
atom at 5 Å into the barrier will feel the same force as at 1.0 Å.
denotes the temperature of the cavity walls in Kelvin. If the system
quasi-temperature is below this setpoint, particles will be accelerated on their
return to the interior. Alternately, they will be retarded if the system is too
warm. A temperature of 0.0 K will turn off wall temperature control, returning
molecules to the system with the same momentum as when they encountered
the barrier.
specifies and/or automatically generates atomic distance constraints. The optional flag,
angstroms, can be used to indicate that data will be entered in Ångstrøms rather than
Bohr.
is the convergence criterion for application of constraints. All distances must
be within +/- tolerance of the specified constraint. Additionally, the RMS
deviation of all constrained distances must be below 2/3 of tolerance.
is the fraction of the total distance correction to be applied on each constraint
iteration.
commands frog to find the closest A atom to each atom X that is closer than
rmax and apply const. The first type line above examines each H atom and
looks for any O atoms within 1.2 Å. The shortest distance, if any, is then fixed at
0.9 Å. Similarly, the second type line binds each F to the closest C within 1.7 Å,
but since const=0.0, that distance is fixed at the current value. The third type
line attaches H atoms to the appropriate nearby C, but at the current average
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:
In this example, starting from the time 2000.0 a.u., velocities are to be scaled every step to keep average total energy constant. Then, from 2500.0 a.u., gradual cooling at the default rate (annealing) is to occur until the time 3000.0 a.u., when free Newtonian dynamics will resume.
Here are all the possible instructions:
These commands cause velocities to be scaled so as to keep the average kinetic energy (i.e. 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).
At some time during the ab initio MD run the user can specify a new value for one of the dynamical variables. The old value is discarded. Single values are given by x=real number. Vectors must be read in frog format from file=file.
In Simulated Annealing MD, the temperature of a run is lowered so as to find 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).
Finally, this instruction turns off any previous action and resumes free dynamics. This is the default status of an MD run.
This keyword allows to carry out Tully-type fewest switches surface hopping (SH) [238]. This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [239]; for the current implementation see Tapavicza et al. 2011 [240]. In the current implementation the surface hopping algorithm only allows switches between the first excited singlet state and the ground state. However, total energies of higher excited states can be computed during the MD simulation. The proper functioning of SH has only been tested for the option
To carry out SH dynamics simulations, the keyword $surface_hopping has to be added to the control and mdmaster file. In addition several keywords are required in the control file:
needed to compute non-adiabatic couplings; this keyword requires the use of
weight derivatives in section dft
keyword needed to collect Cartesian non-adiabatic coupling vectors along the
trajectory
keyword needed to ensure dynamics starting in S1
collects time integration of excitation energies along the trajectory
collects amplitudes of the adiabatic states along the trajectory
Special caution has to be taken to control problems related to conical intersections [241,242]. At geometries where conical intersections between the ground and excited state are present, DFT often exhibits singlet instabilities, which leads to imaginary excitation energies in linear response TDDFT; in this case the MD run is terminated. This problem can be circumvented by the use of the Tamm-Dancoff approximation (TDA) to TDDFT (see 8). In addition an optional keyword for the md_master file can be used:
enforces a switch to the ground state in case the 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.
In order to control the program execution, you can use the following keywords within the control file:
Switches on the calculation of the MP2 NMR shieldings. The required SCF shielding
step will be performed in the same run. This flag will be set by the script mp2prep.
specifies the number of loops (or ’passes’) over occupied orbitals n when doing an
MP2 calculation: the more passes the smaller file space requirements—but CPU time
will go up. This flag will be set by the script mp2prep.
Scratch file settings for an MP2 calculation. Please refer to Section 20.2.18 for a
description of the syntax. This flag will be set by the script mp2prep.
Sets the convergence threshold for the shielding constant of succeeding CPHF
iterations. The unit is ppm and the default value is 0.01.
This selects the atom number for convergence check after each cphf iteration. After
this convergence is reached all other atoms are checked, too (default: 1).
have the save meaning as in dscf (see Section 20.2.7);
Since mpshift works ’semi-direct’ it uses the same integral storage.
The scratch files allocated by mpshift can be placed anywhere in your file systems instead
of the working directory by referencing their pathnames in this data group. All possible
scratch files are listed in the following example:
stands for traloop start and traloop end. Each loop or pass in MP2 chemical shift
calculations can be done individually by providing the keywords $trast and $trand. This
can be used to do a simple parallelization of the run:
Create separate inputs for each traloop. Add
in the control files, number goes from 1 to the number of $traloops. Each calculation will create a restart file called restart.mpshift. To collect all steps and to do the remaining work, copy all restart files to one directory and rename them to restart.mpshift.number, add $trast -1 and $trand number_of_traloops to the control file and start mpshift.
UaiBβ will be written to file for a subsequent calculation of VCD rotational strengths by
aoforce.
Note: All options in this section are obsolete since the task distribution is done automatically — but they are kept for compatibility reasons.
On all systems the parallel input preparation is done automatically. Details for the parallel installation are given in Section 3.2.1. The following keywords can be used for parallel runs, but are all optional:
$parallel_platform architecture
Currently the following parallel platforms are supported:
for systems with very fast communication; all CPUs are used for the linear
algebra part. Synonyms for SMP are:
HP V-Class, SP3-SMP and HP S/X-Class
for systems with fast communication like Fast-Ethernet, the number of CPUs
that will be taken for linear algebra part depends on the size of the matrices.
Synonyms for MPP are:
SP3 and linuxcluster
for systems with slow communication, the linear algebra part will be done on
one single node. Synonyms for cluster are:
HP Cluster and every platform that is not known by TURBOMOLE
If you want to run mpgrad, $traloop has to be equal to or a multiple of the number of parallel workers.
For very large parallel runs it may be impossible to allocate the scratch files in the working directory. In this case the $scratch files option can be specified; an example for a dscf run is given below. The scratch directory must be accessible from all nodes.
For DFT with dscf/grad$pardft can be specified:
The tasksize is the approximate number of points in one DFT task (default: 1000) and memdiv says whether the nodes are dedicated exclusively to your job (memdiv=1) or not (default: memdiv=0).
For dscf and grad runs, a parallel statistics file which has to be generated in advance can be used.
The filename is specified with
$2e-ints_shell_statistics file=DSCF-par-stat
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:
and for a grad statistics run:
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.