25.2.10 Keywords for Modules `dscf` and `ridft`

$denconv real

Convergency criterion for the root mean square of the density matrix. If you want to calculate an analytical MP2 gradient (program mpgrad) real should be 1.d-7 or less. By default, a criterion of 1.d-6 is applied for Hartee–Fock calculations.

$dft options

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

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

functional name

Specification of the functional, default is BP86, printed as functional b-p. For all possible—and useful—functionals, please refer to page 25.2.10 and for definition of the functionals the section 6.3 on page 6.3.

Example (default input):

$dft
   functional b-p

gridsize integer or minteger

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

Example:

$dft
   gridsize m3

gridtype integer

—not recommended for use—
Specification of the mapping of the radial grid.

Possible values for gridtype are $1, \ldots, 6$, but gridtype 4 to 6 is only for the use with functional lhf (see page 25.2.10). For the definition of gridtype 1 – 3, please refer to Eq. (16), (17) and, (19) in Ref. [416].

Example (default value):

$dft
   gridtype 3

debug integer

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

nkk integer

Specification of the sharpness of the partition function as proposed by Becke [417], default is nkk 3. The usage of nkk makes sense only in the range $1 \le$ nkk $\le 6$.

Example:

$dft
   nkk 3

—not recommended for use—
Only for user-specified Lobatto grids, i.e. gridsize 9, ntheta specifies the number of $\theta$ points and nphi specifies the number of $\phi$ points. For the fixed Lobatto grid, i.e. gridtype 8, the default value is ntheta 25 and nphi 48.

When gridsize 9 is given, you have to specify both, ntheta and nphi (see below), otherwise the program will crash. The restriction for user-defined Lobatto grids is the number of grid points, which must not exceed 2000 grid points.

Example:

$dft
  gridsize 9
  ntheta 30
  nphi 60

old_RbCs_xi

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

Example:

$dft
   old_RbCs_xi

radsize integer

Specifies the number of radial grid points. Default values depend on type of atom and grid (see keyword gridsize). The formula for the radial gridsize is given as, \[\begin{equation*} \mathrm{number\ of\ radial\ grid\ points} = \mathrm{ioffrad} + (\texttt{radsize}-1)*5\,. \end{equation*}\]

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 grids for relativistic calculations, i.e. gridsize 3a, use the atom number Z in a modified formula given as, \[\begin{equation*} \mathrm{number\ of\ radial\ grid\ points} = 20 + (\texttt{radsize}-1)*5 + Z\,. \end{equation*}\]

The radial grid size increases further for finer grids:

`gridsize`	1	2	3	4	5	6	7	8	9
`radsize`	1	2	3	6	8	10	14	9	3

If you want to converge results with respect to radial grid size, increase radsize in steps of 5, which is convenient (see equation above).

diffuse integer

Serves to increase quadrature grids; this is recommended in case of very diffuse wavefunctions. With the keyword diffuse grids are modified by changing the linear scaling factor $\xi$ of radial grid points and the radial gridsize:
radsize $\Longrightarrow$ radsize + incr
$\xi$ $\Longrightarrow$ $\xi$ * 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 $\xi$, see Eq. (16)–(19) and Table 1 in Ref. [416].

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

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

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

Example (Recommendation):

$dft
  gridsize m4
  diffuse 2

—for developers only—
Radial grid points have a linear scaling parameter $\xi$, see Eq.(16)–(19) and Table 1 in Ref. [416]. With the following input,

$dft
   rhostart 50
   rhostop 200

one performs a numerical integration for the density and the exchange correlation term for $\xi=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 $\xi$ values, see Figure 1 in Ref. [416]. Program stops after this testing.

reference

Usage of the reference grid, which is a very fine grid with very tight thresholds. The default values for the different variables are:

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

Please refer to the corresponding sub-keywords for explanation.

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

$dft
   functional b-p
   reference

Note: The sub-keywords radsize, fcut, dgrenze, fgrenze, qgrenze can be used to overwrite the above values, but not vice versa.

test-integ

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

batchsize integer

Grid points are sorted into batches, which are then processed. This increases efficency. This should be changed only by developers. Default is batchsize 100.

fullshell

Standard grids have reduced number of spherical grid points near nuclei. With the keyword fullshell this reduction is suppressed. Reference grid (see keyword reference) always has full spherical grids with 1202 points. Should be used to checked the influence of spherical grid reduction.

Example for the usage of fullshell:

$dft
   functional b-p
   gridsize m4
   fullshell

—for developers only—
Values of real effects efficiency of the quadrature, default is symblock1 0.001 and symblock2 0.001, one can try higher or smaller values.

xparameter integer

—not recommended for use—
Where xparameter (default) can be: sgrenze (8), fgrenze (10), qgrenze (12), dgrenze (12) and, fcut (14). These parameters control neglect of near zeros of various quantities. With xparameter integer one changes the default. integer larger than defaults will increase the numerical accuracy. Tighter threshold are set automatically with keyword $scfconv (see section 25.2.10 on page 25.2.10).

weight derivatives

Includes the derivatives of quadrature weights to get more accurate results. Default is that the derivatives of quadrature weights will be not considered, see section 25.2.16 on page 25.2.16.

gridordering

Grid points are ordered into batches of neighbouring points. This increases efficiency, since now zeros can be skipped for entire batches. gridordering is default for serial version, not for the parallel one. You cannot use weight derivatives and gridordering together.

Example for switching off gridordering:

$dft
   gridordering 0

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 the neglection of certain shell pairs and thus accelerates the calculation. s-junc and p-junc set the corresponding thresholds to $\mathrm{10}^{-\mathit{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 for p-junc 6). The defaults settings are s-junc 6 and p-junc 6 for grad, rdgrad, and ridft. Too small values for s-junc and p-junc might result in non-converging calculations or deteriorated results. This holds particularly for excited state calculations. For escf the default values are therefore s-junc 8 and p-junc 8.

$electrostatic field

Specification of external electrostatic field(s). The specification may take place either by Ex, Ey, Ez or by x, y, z, |E|. See also $fldopt.

Example:

$electrostatic field
       0.1000E-03   0.000       0.000

Electrostatic fields are also possible for 2c calculations, turning on geofield is mandatory in these calculations. $pcc activates picture change correction for fields integrals in 1c and 2c (local) BSS/DKH/X2C calculations.

$magnetic field

Specification of a single external magnetic field. The specification may take place either by Bx, By, Bz or by x, y, z, |B|.

Example:

$magnetic field
       0.000   0.000   0.1000E-03  
       tesla
       z-field
       nomag

Magnetic fields are currently only possible for 2c calculations.

The optional keyword $tesla specifies the field is given in Tesla. Otherwise, atomic units are assumed.

The optional keyword $z-field specifies that only the magnetization in z-direction is significant (i.e. no broken symmetry solutions).

The optional keyword $nomag specifies that the magnetization densities vanish. This is usually true for small fields and closed shell references. This option will enforce TD-DFT to be possible also in magnetic fields.

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

Requests calculation of occupation numbers at a finite temperature $T$. For an orbital with the energy $\varepsilon _{i}$ the occupation number $n_{i}\in \left[ 0,1\right]$ is calculated as \[n_{i}=\frac{1}{2}\text{erfc}\left( \frac{\varepsilon _{i}-\mu }{fT}\right) \text{,}\] where $\mu$ is the Fermi level. The factor $f=4k/\sqrt{\pi}$ is chosen to yield the same slope at the Fermi level as the Fermi distribution [55].

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

$firstorder

Perform first-order SCF-calculation, i.e. perform only one SCF-iteration with the start MOs (which should be the orthogonalized MOs of two independent subsystems as is explained in detail in Chapter 22).

$fldopt options

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

1st derivative on/off: Calculate numerical 1st derivative of SCF energy with respect to electrostatic field (default: off), increment for numerical differentiation is edelt (see below).
2nd derivative on/off: Calculate numerical 2nd derivative of SCF energy with respect to electrostatic field (default: off), increment for numerical differentiation is edelt.
edelt= real: Increment for numerical differentiation (default: 0.005).
fields on/off: Calculate SCF energy for non-zero external electrostatic
fields defined in $electrostatic field.
geofield on/off: Calculate SCF energy for one external field definition and dump disturbed MOs onto $scfmo. This enables to evaluate properties or perform geometry optimizations in the presence of an external field.

Caution: don’t use the RI approximation for all these calculations since this will lead to non-negligible errors!!

$hsout

Invokes writeout of Fock matrix and overlap matrix for a converged ridft calculation. For open-shell calculations, two Fock matrices are written.

$incore integer

By using this option the two-electron integrals are kept in RAM; integer specifies how many megabytes should be allocated. If the integrals exceed the RAM allocated the program reverts to the standard mode. Supports all methods which process two-electron integrals, i.e. SCF and DFT (including hybrid functionals); RHF and UHF.

The following condition must be met:

$scfdenapproxl 1

and rhfshells 1 or 2. It is advisable to set $thize as small as possible (e.g. $thize 0.1d-08) and to remove the keyword $scfdump.
Note: this keyword does not work for parallel runs.

$intsdebug sao

Output of one-electron matrices expressed in symmetry-adapted AO basis. Can be used in C1 symmetry or with symmetry. In case of symmetry, the matrices are printed for each irreducible representation. Note that $intsdebug cao is not supported anymore.

Note that the output gives one triangle of the one-electron matrices. Thus the entries are:

(11) (21) (22)
(31) (32) (33)
(41) (42) (43)
…

In C1, the order of the basis functions is such that all s-orbitals are given first, then all p-orbitals, all d-orbitals and so on. So we have:

1. atom 1s,2s,3s…
2. atom 1s,2s,3s…
…
Note that the SAO basis uses 5 d-functions and CAO uses 6.

$mo-diagram only nirreps=integer

If this keyword is set the energies and symmetry labels of all occupied MOs will be dumped to this data group. This may be helpful to draw mo-diagrams. If only has been set only the start MOs are dumped and the program quits.

nirreps will hold the total number of displayed orbitals after the successful run.

$mom

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

$moprint

If this keyword is present all occupied orbitals are dumped to standard output. Be careful about this option as it can create huge output files in case of many basis functions.

$mo output format format

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

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

$natural orbitals

This data group will be written after an UHF calculation (together with $natural orbital occupation) and contains the natural space orbitals (same syntax as $scfmo).

$natural orbital occupation

This data group will be written after an UHF calculation (together with $natural orbitals) and contains the occupation of natural orbitals (syntax as any data group related with orbital occupation information, e.g. $closed shells), e.g.

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

$prediag

concerns the first SCF iteration cycle if start MOs from an EHT guess are used.

The SCF iteration procedure requires control mechanisms to ensure (fast) convergence, in TURBOMOLE these are based on orbital energies $\epsilon_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 $\epsilon_i$ are not available. The keyword $prediag provides ’$\epsilon_i$ of the zeroth iteration’ by diagonalization of occ–occ and virt–virt part of the first Fock matrix, to allow for level shifting etc.. See $scfdiis below.

$restart dscf twoint

Try a dscf restart. The program will read the data group $restartd (which must exist, also $scfmo has to exist!) and continue the calculation at the point where it ended before. If the additional option twoint is appended, the program will read the two-electron integrals from the files specified in $scfintunit, so there will be almost no loss of cpu-time.

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

$restartd data

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

$rundimensions

is set by define so usually no changes are necessary. The dimensions must be greater or equal to those actually required, i.e. you can delete basis functions and keep rundimensions. This keyword is not necessary for small cases.
Example:

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

$scfconv integer

SCF convergency criterion will be $10^{\mathit{-integer}}$ for the energy. Gradients will only be evaluated if integer $>$ 6.

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

Damping parameters for SCF iterations in order to reduce oscillations. The old Fock-operator is added to the current one with weight 0.5 as start; if convergence is good, this weight is then reduced by the step 0.05 in each successive iteration until the minimum of 0.1 is reached. (These are the default settings of define for closed-shell RHF). DSCF automatically tries to adjust the weight to optimize convergence but in difficult cases it is recommended to start with a large weight, e.g. 1.5, and to set the minimum to a larger value, e.g. 0.5.

$scfdebug options

Flags for debugging purposes. Following options are available:

vectors integer

Output level concerning molecular orbitals. integer=0 (default) means minimal output,

$>$1 will output all start MOs and all MOs in each iteration.

density integer

Output level concerning difference density matrices.

debug integer

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

$scfdenapproxl integer

Direct SCF procedures build the Fock matrix by exploiting information from previous iterations for better efficiency. With this keyword information from the last integer iterations will be used. This feature is switched on with the default value 20, even if the keyword is absent. The user may reduce the number of iterations employed to smaller values (e. g. 10) in cases were numerical stability could become an issue. With the value 0 this feature is switched off; the Fock matrix is constructed from scratch in each iteration.

$scfdiis options

Control block for convergence acceleration via Pulay’s DIIS ¹.
Options are:

errvec=char

specifies the kind of error vector to be used (two different kind of DIIS algorithms)

char=’FDS’ or ’SDF’ or ’FDS-SDF’: uses $\rm (FDS-SDF)$ as error vector.
char= none: no DIIS
char= sFDs: use $S^{-1/2}FDS^{1/2}-\mathrm{transposed}$

Further suboptions:

maxiter=integer

maximum number of iterations used for extrapolation.

debug=integer

debug level (default: 0)

integer=1: print applied DIIS coefficients
integer=2: print DIIS matrix and eigenvalues, too

qscal=real

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

thrd=real

directs the reduction of qscal to qscal = 1.0 (no damping in DIIS), done if $||errvec||$ $<$ thrd.

Defaults for $prediag (see above) and $scfdiis

errvec=FDS-SDF, maxiter=5, qscal=1.2, thrd=0.0, this implies DIIS damping in all iterations, prediag is switched of.

Recommended:

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

$scfdump

Dump SCF restart information onto data group $restartd and dump SCF MOs in each iteration onto $scfmo (scfdump = iter). Additionally, a data block $scfiterinfo will be dumped containing accumulated SCF total-, one- and two-electron energies of all previous SCF iterations. Information that will allow you to perform a restart if your calculation aborts will be dumped on data group $restartd (see also $restart).

$scfintunit options

Disc space specification for two-electron integrals. The following suboptions are available (and necessary):

unit=integer: A legacy suboption which is ignored.
size=integer: Filespace in megabytes for this file. size=0 leads to a fully direct run. size is set by a statistics run, see $statistics. DSCF switches to direct mode if the file space is exhausted.
file=char: Filename. This may also be a complete path name, if you want to store the integrals in a special directory. Make sure the file is local, otherwise integrals are transmitted over the network.

Thus your data group $scfintunit may look like this:

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

Maximal 30 files may be specified in this way.

$scfiterlimit integer

Maximum number of SCF iterations (default: 30).

$scfmo none file=char

Input/output data group for SCF MOs. You can specify

none: To perform a calculation without a start vector (i.e. use a core Hamiltonian guess).
file=char: The file where the MOs are written on output (default: mos).

These two options can also be used for $uhfmo_alpha and $uhfmo_beta to use a core guess and write the molecular orbitals to file.

After running define or a TURBOMOLE calculation additional options may appear specifying the origin of the MOs:

expanded: These MOs were obtained by projection form another basis set. They should not be used for wavefunction analysis.
scfconv=integer: The MOs are converged SCF MOs , the convergence criterion applied was $10^{\mathit{-integer}}$
scfdump=integer: The MOs are unconverged SCF MOs which were written on this data group after iteration integer. The latter three options are mutually exclusive.
format(format string): This specifies the FORTRAN format specification which was used for MO output. The standard format is (4d20.14). (See data group $mo output format.)

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

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

$scforbitalorder on/off

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

$scforbitalshift options

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

Options are:

noautomatic

Automatic virtual shell shift switched off.

automatic real

Automatic virtual shell shift switched on; the energies of virtual orbitals will be shifted if the HOMO-LUMO gap drops below real such that a gap of real is sustained. This is the default setting if the keyword is missing with real=0.1.

closedshell=real

Option for open-shell cases. Closed shells are shifted to lower energies by real. The default shift value is closedshell=0.4.
Note: Normally this will disable the automatic shift of energies of virtual orbitals! To override this, you should append an exclamation mark to the ’automatic’ switch, i.e. specify ’automatic! real’.

individual

Set shifts for special occupied MOs. To use this option, start the line with the symmetry label and the list of MOs within this symmetry and append the desired shift in brackets as in the following example:

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

$scftol real

Integral evaluation threshold. Integrals smaller than real will not be evaluated. Note that this threshold may affect accuracy and the convergence properties if it is chosen too large. If $scftol is absent, a default value will be taken obtained from $scfconv by $real=\frac{10^{-(scfconv + 1)}}{3 \cdot \#bf}$ (#bf = number of basis functions).

$scratch files

The scratch files allocated by dscf can be placed anywhere in your file systems instead of the working directory by referencing their path names in this data group. All possible scratch files are listed in the following example:

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

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

$statistics options

The following options are allowed

off: Do not perform integrals statistics
dscf: Perform integrals statistics for dscf
mpgrad: see mpgrad
dscf parallel: see PARALLEL PROCESSING

Options dscf parallel, grad, mpgrad will be described in the related chapters.
If $statistics dscf has been given integral prescreening will be performed (which is an $n^4$-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 files pace needed for the current combination of $thize and $thime will be written to the data group ($scfintunit) and $statistics dscf will be replaced by $statistics off.

$thime

Integral storage parameter, which is related to the time needed to calculate the integral. The larger integer the less integrals will be stored. The default value is integer = 5. (see also $thize, $statistics)

$thize

Integral storage parameter, that determines, together with $thime, the number of integrals stored on disc. Only integrals larger than real will be stored. The default value is real = 0.100E-04.

RHF/ROHF

$closed shells

Specification of MO occupation for RHF, e.g.

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

$open shells type=1

MO occupation of open shells and number of open shells. ’type=1’ here means that there is only a single open shell consisting e.g. of two MOs:

b2g     1                                      ( 1 )
b3g     1                                      ( 1 )

$rohf

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

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

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

$roothaan

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

$roothaan
    a = 3/4      b = 3/2

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

For further information on ROHF calculations (e.g. with more than one open shell), see the sample input in Section 26.6 and the tables of Roothaan parameters in Section 6.4.

Note that this keyword toggles the ROHF mode also for more than one open shell. If it is not given, the open-shell electrons are simply ignored.

UHF

$alpha shells and $beta shells

these two data groups specify the occupation of alpha and beta spin UHF MOs (syntax as any data group related with orbital occupation information, e.g. $closed shells)
Example:

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

$uhf

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

$spin constraint {real}

The presence of $uhf activates the UHF mode. An additional constraint can be imposed for the (spin)$^2$ exceptation value by using $spin constraint {real}. A non-zero real value fixes the (spin)$^2$ expectation value by using [84]
$\rm F^{\alpha}(SAO) - 2*\tau*SD^{\beta}S$
$\rm F^{\beta}(SAO) - 2*\tau*SD^{\alpha}S$
instead of $\rm F^{\alpha},F^{\beta};$ in the limiting case (real $\rightarrow$ INFINITY) the ROHF (spin)$^2$ expectation value would result. Default spin constraint 0, this is also used if the key is not set. We suggest to increase the constraint in the sequence 0.01, 0.1, 1.0, and 10.0 to check the impact on the (spin)$^2$ expectation value.

$uhfmo_alpha and $uhfmo_beta

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

$uhfmo_beta

see $uhfmo_alpha

DFT

$dft
 functional b-p
 gridsize m3

for DFT calculations one has to specify the functional and the grid (for the quadrature of the exchange correlation part). The settings above are default: both lines can be left out if the B-P86 functional and grid m3 are required. Other useful functionals supported are:

b-lyp
b3-lyp
b3-lyp_Gaussian: (equivalent to the Gaussian98 keyword B3LYP with VWNIII)
bh-lyp
s-vwn
s-vwn_Gaussian: (equivalent to the Gaussian98 keyword SVWN with VWNIII)
tpss
tpssh

Possible grids are 1–7 and m3–m5 where grid 1 is coarse (least accurate) and 5 most dense. The so-called multiple grids m3–m5 use the following ansatz: 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 general response properties and parallel calculations, the multiple grids are not recommended, see also Sec. 8. Multiple grids should not be used for magnetic properties. In relativistic all-electron calculations grids with an increased number of radial points are recommended. These can be selected by appending the character a after the gridsize, e.g. gridsize 4a. For a reference calculation with a very fine grid and very tight thresholds use ’reference’ as grid specification instead of ’gridsize xy’.

Note: the functionals b3-lyp_Gaussian and s-vwn_Gaussian are made available only for comparability with Gaussian. The functional VWNIII is much less well founded than VWN5 and the TURBOMOLE team does not recommend the use of VWNIII.

RI

Dscf does not run with the keyword $rij: you must call the RI modules Ridft and Rdgrad for energy and gradient calculations. However, it does run with the keyword $rik, but it will ignore all RI settings and do a conventional non-RI Hartree–Fock or DFT calculation.

$rij: Enforces an RI-$J$ calculation if module ridft is used, can be used for Hartree-Fock as well as for DFT calculations with pure or hybrid functionals.
$ridft: Obsolete keyword - use $rij instead!
$rik: Enforces a RI-JK calculation if module ridft is used, can be used for Hartree-Fock as well as for DFT calculations with pure or hybrid functionals.
$ricore integer: Choose the memory core available (in megabyte) for special arrays in the RI calculation (the less memory you give the more integrals are treated directly, i.e. recomputed on the fly in every iteration)
$jbas file=auxbasis: Cross reference for the file specifying the auxiliary basis as referenced in $atoms. We strongly recommend using auxbasis sets optimized for the respective MO basis sets, e.g. use SVP (or TZVP) for the basis and the corresponding auxbasis as provided by define (default: file=auxbasis).
$ripop: Calculation of atomic charges according to the $s$ partial wave and atomic dipole moments according to the $p$ partial wave as resulting from the auxbasis representation of the density

RI-JK

If the keyword $rik is found in the control file, ridft performs a Hartree–Fock–SCF calculation using the RI-approximation for both Coulomb and HF-exchange (efficient for large basis sets). For this purpose needed (apart from $ricore):

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

MARI-J

Multipole-Accelerated-Resolution-of-Identity-$J$. This method partitions the Coulomb interactions in the near- and far-field parts. The calculation of the far-field part is performed by application of the multipole expansions and the near-field part is evaluated employing the RI-$J$ approximation. It speeds up calculation of the Coulomb term for large systems. It can only be used with the ridft module and requires setting of the $ridft keyword.

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

The following options are available:

precision: specifies precision parameter for the multipole expansions. Low-precision MARI-J calculations require $1\cdot 10^{-6}$, which is the default. For higher precision calculations it should be set to $1\cdot 10^{-8}$–$1\cdot 10^{-9}$.
lmaxmom: maximum l-moment of multipole expansions. It should be set to a value equal at least twice the maximum angular momentum of basis functions. Default value is 10 and it should probably never be set higher than 18.
thrmom: Threshold for moment summation. For highly accurate calculations it should be set to $1\cdot 10^{-24}$.
nbinmax: number of bins per atom for partitioning of electron densities. Default value is 8 and hardly ever needs to be changed.
wsindex: minimum separation between bins. Only bins separated more than the sum of their extents plus wsindex are considered as far-field. Default is 0.0 and should be changed only by the experts.
extmax: maximum extent for charge distributions of partitioned densities. Extents with values larger then this are set to extmax. Hardly ever needs to be changed.

Seminumeric HF-Exchange

If the keyword $senex is found in the control file, ridft performs a Hartree–Fock– SCF calculation using the seminumerical approximation for HF-exchange. Standard dft-grids can be used for the numerical integration. Smaller grids (-1,0) and the corresponding m-grids (m1,m2) have been defined as well. For high precision energies the use of de-aliasing is recommended (do_sfit, C1 symmetry only) which will yield reliable energies with the m1 grid in nearly all cases. Alternatively grids of at least size m3 are recommended for heavy atoms. The gridsize can be modified just like in dft-calculations. We introduced three different keywords to finetune where you want to employ the seminumerical exchange approximations. The keyword $senex activates seminumerical calculations in energies (ridft), excitation energies (escf), vibrational frequencies (aoforce) and chemical shifts (mpshift). The keyword $esenex activates seminumerical exchange in escf, aoforce, mpshift, egrad(excitation energy part only) but not in ridft. Further it resets the default grid to -1 in these modules; which was found to be sufficient in most cases. The keyword $dsenex activates seminumerical gradient calculations in ground (rdgrad) and excited states (egrad; excitation energy AND gradient part). An example using the default grid and de-aliasing for SCF (m1) and grid m2 for gradients looks like this:

$senex
 do_sfit
$dsenex
 gridsize m2

Seminumeric Coulomb+Exchange: The pseudospectral approach

The modules escf, aoforceand egradcan make use of the pseudospectral approximation to calculate Coulomb contribution seminumerically on a grid. To do so simply add the keyword $pseudospectral additionally if $senex or $esenex are present. This is especially valuable when no RI-Fitting basis is available for the chosen basis set (e.g. "cbas" for Sapporo-type basis sets, ANO-type basis sets, x2c-TZVPPall etc.) as the pseudospectral approach does not need any auxiliar/fitting basis sets at all. Excitation energies, harmonic frequencies and ZPEs obtained from the pseudospectral approximation are extremely accurate (usually considerably lower errors than RI) and still orders of magnitude faster than using classic Coulomb integrals. For ground states (e.g. ridft) the pseudospectral approximation is not yet available. As the RI approximation provides true upper bounds in this case it will outperform the pseudospectral approximation. Further auxiliary basis sets of type "jbas" are, unlike "cbas" auxiliary basis sets, widely available for the whole periodic table. In C1-symmetry we recommend the use of de-aliasing using the do_sfit keyword. An example for applying the de-aliased pseudospectral approximation looks like this:

$senex
 do_sfit
$pseudospectral

LHF

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

$dft
 functional lhf
 gridsize 6

With the LHF potential Rydberg series of virtual orbitals can be obtained. To that end, diffuse orbital basis sets have to be used and special grids are required.

gridtype 4 is the most diffuse with special radial scaling; gridtype 5 is for very good Rydberg orbitals; gridtype 6 (default in Lhfprep) is the least diffuse, only for the first Rydberg orbitals.

Only gridsize 3–5 can be used, no multiple grids.

Use test-integ to check if the selected grid is accurate enough for the employed basis-set.

How to do LHF runs

Do a Hartree–Fock calculation using dscf.
Use the script lhfprep to prepare the control file (the old control file will be saved in control.hf and the molecular orbitals in mos.hf or in alpha.hf and beta.hf for the spin-unrestricted case). See lhfprep -help for options.
Run again dscf.

Otherwise the LHF functional can be selected in define: in this case default options are used.

Options for the LHF potential can be specified as follows ( see also lhfprep -help)

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

off-diag off

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

numerical-slater on

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

numerical-slater off

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

asymptotic

for asymptotic treatment there are three options:

asymptotic off: No asymptotic-treatment and no use of the numerical Slater. The total exchange potential is just replaced by $-1/r$ in the asymptotic region. This method is the fastest one but can be used only for the density-matrix convergence or if Rydberg virtual orbitals are of no interest.
asymptotic on: Full asymptotic-treatment and use of the numerical Slater in the near asymptotic-region.
asymptotic dynamic=1.d-3: Automatic switching on (off) to the special asymptotic treatment if the differential density-matrix rms is below (above) 1.d-3. This is the default.

pot-file save

the converged Slater and correction potentials for all grid points are saved in the files slater.pot and corrct.pot, respectively. Using pot-file load, the Slater potential is not calculated but read from slater.pot (the correction potential is instead recalculated). For spin unrestricted calculations the corresponding files are slaterA.pot, slaterB.pot, corrctA.pot and correctB.pot.

homo

allows the user to specify which occupied orbital will not be included in the calculation of correction potential: by default the highest occupied orbital is selected. This option is useful for those systems where the HOMO of the starting orbitals (e.g. EHT, HF) is different from the final LHF HOMO. homob is for the beta spin.

correlation func=

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 $R_F \; \; \Delta_F$
$0 \ldots R_F-\Delta_F$ : basis-set correction potential
$R_F-\Delta_F \ldots R_F+\Delta_F$ : smooth region
$R_F+\Delta_F \ldots +\infty$ : asymptotic correction
Defaults: $R_F =10$, $\Delta_F=0.5$
: slater-region $R_N \; \; \Delta_N \; \; R'_F \; \;\Delta'_F$
$0 \ldots R_N-\Delta_N$ : basis-set Slater potential
$R_N-\Delta_N \ldots R_N+\Delta_N$ : smoothing region
$R_N+\Delta_N \ldots R'_F-\Delta'_F$ : numerical Slater
$R'_F-\Delta'_F \ldots R'_F+\Delta'_F$ : smoothing region
$R'_F+\Delta'_F \ldots +\infty$ : asymptotic Slater
Note: $R'_F-\Delta'_F \le R_F-\Delta_F$
Defaults: $R_N=7$, $\Delta_N=0.5$, $R'_F =10$, $\Delta'_F=0.5$
Use correct-b-region and slater-b-region for the beta spin.

Two-component SCF (GHF)

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

$soghf: enforces two-component-SCF calculations; this option is compatible with $rij, $rik and $dft.
$kramers: switches on Kramers-restricted formalism
$collinear: switches on collinear two-component formalism (not rotational invariant)
$gdiis: enforces DIIS for complex Fock operator.

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

Relativistic all-electron calculations can be done employing the X2C, the BSS or the DKH Hamiltonian. Implemented for modules dscf and ridft. Note that gradients are only available with X2C in the modules grad and rdgrad.

$rx2c

switches on X2C calculation.

$rbss

switches on BSS calculation.

$rdkh integer

switches on DKH calculation of order integer.

$dkhparam integer

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

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

Note in particular that the parametrization does not affect the Hamiltonian up to 4th 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 [419].

$rlocal

switches on local approach. Default is DLU (see $rlocpara).

$rlocpara integer

selects parameterization of the local approximation. Valid values are 0, 1, or 2. Default is the DLU (0). 1 refers to the DLH approximation. 2 applies the low-scaling DLU variant, which approximates the atomic-off diagonal block of distant atoms with the non-relativistic approximation of the respective balance condition. For details on the different parametrizations see [164]. Gradients are restricted to the DLU approximations.

$finnuc

switches on finite nucleus model based on a Gaussian charge distribution with parameters taken from [353]. In case of magnetic properties, the vector potential is also based on finite nuclei.

$snso

selects screened-nuclear-spin-orbit (SNSO) approach for the spin-dependent integrals in two-component calculations.

$snsopara integer

selects parametrization of SNSO. Valid values are 0-6. 0 refers to the original set of parameters originally used in low-order DKH theory [167], while 1 refers to the modified parameters [168, 169] for exact decoupling methods. Default is 1. 2 and 3 are taken from [340] and [341]. These can only be used with the Pauli Hamiltonian and the non-relativistic framework. Options 4 and 5 are taken from [190] and describe the Dirac–Coulomb and the Dirac–Coulomb–Breit approach. 6 applies the later parameters of [191].

$pcc

switches on relativistic picture-change correction for the calculation of expectation values.

All of these keywords are compatible with $soghf employing the RI-$J$ approximation.

P.Pulay; Chem. Phys. Lett., 73, 393 (1980), P.Pulay; J. Comput. Chem., 4, 556 (1982)↩︎
C. C. J. Roothaan: Rev. Mod. Phys. 32 (1960) 179.↩︎

25.2.10 Keywords for Modules dscf and ridft

RHF/ROHF

UHF

DFT

RI

RI-JK

MARI-J

Seminumeric HF-Exchange

Seminumeric Coulomb+Exchange: The pseudospectral approach

LHF

How to do LHF runs

For expert users

Two-component SCF (GHF)

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

25.2.10 Keywords for Modules `dscf` and `ridft`