8.4 How to Perform

The most convenient way to set up an escf or egrad calculation is to use the ex option of the last (“general”) define menu, see Chapter 4. define will automatically provide most of the keywords discussed below.

A large number of (not necessarily realistic) sample inputs is contained in the escf and egrad subdirectories of the test suite (TURBOTEST directory).

8.4.1 Preliminaries

All response calculations require a complete set of converged (occupied and virtual) SCF MOs. It is strongly recommended to use well converged MOs, since the error in the ground-state wavefunction enters linearly in all response properties. Thus, before starting escf or egrad, specify the keywords

$scfconv 7
$denconv 1d-7

in control, perform a dscf statistics run if semi-direct integral processing is to be used (see Chapter 3), and (re-)run dscf or ridft,

dscf > dscf.out &

or

ridft > ridft.out &

in case of RI-J.

The above tight convergence criteria are also recommended for excited state geometry optimizations. It is also recommended to avoid multiple grids as they negatively influence the numerical stability (use 3 instead of m3). To perform a two-component TDDFT calculation, the two-component version of ridft has to be run before (see Chapter 6.5) using the keywords $soghf and $kramers.

8.4.2 Polarizabilities and Optical Rotations

The calculation of dynamic polarizabilities is controlled by the keyword

unit specifies the unit of the following frequencies and may be ev, nm, 1/cm, or a.u. (default). The frequencies may be either purely real or purely imaginary. For example, to calculate dynamic polarizabilities at 590 nm and 400i nm (i is the imaginary unit), specify

$scfinstab dynpol nm
  590
  400 i 

and run escf,

escf > escf.out &

.

The resulting polarizabilities and rotatory dispersions are given in a.u. in the program output (escf.out in the above example).

The conversion of the optical rotation in a.u. to the specific rotation \([\alpha]_{\omega}\) in deg\(\cdot\)[dm\(\cdot\)(g/cc)]\(^{-1}\) is given in Eq. (15) of ref. [223]. \[\begin{equation} [\alpha]_{\omega} = C \cdot \delta(\omega) \end{equation}\](8.26) where \(C = 1.343 \cdot 10^{-4} \omega^2/M\) with \(M\) being the the molar mass in g/mol, \(\omega\) the frequency in cm\(^{-1}\), and \(\delta(\omega)\) is 1/3 trace of the electronic rotatory dispersion tensor given in atomic units.
Please note that \(\delta(\omega)\) has the wrong sign in older TURBOMOLE versions. It has been corrected in version 6.2.

Note that convergence problems may occur if a frequency is close to an electronic excitation energy. This is a consequence of the (physical) fact that the response diverges at the excitation energies, and not a problem of the algorithm.

Static polarizabilities are calculated most efficiently by specifying

$scfinstab polly

before starting escf. This keyword can be combined with $soghf and the other relativistic keywords such as $rdkh, $rbss, $rx2c, $rlocal, and $pcc. Please see Sec. 6.5 for details.

8.4.3 Damped Response Calculations

A damped response calculation (adding a complex constant to a polarizability) can be invoked by additionally adding the following keywords to the control file:

$scfinstab dynpol nm
  590
$damped_response 0.25 eV

The last keyword distinguishes it from a simple polarizability calculation described above. Defined units for the frequency and the damping factor may be different. Valid units are eV, cm-1 and nm; if nothing is specified atomic units are assumed. Further it is possible to specify an evenly spaced list of frequencies:

$scfinstab dynpol eV
  2.0 3.0 10
$damped_response 0.25 eV

This will perform a damped response calculation of 10 frequencies between 2.0\(\,\)eV and 3.0\(\,\)eV with a complex damping factor of 0.25 eV. Damped response polarizabilities can be performed for 1c (open and closed shell) and 2c Kramers-symmetric quasi-relativistic references (ECPs or relativistic all-electron approaches) within the time-dependent DFT (including local hybrid functionals) and also using the Bethe-Salpeter equation. The latter is especially recommended if core states are targeted. Note that in the case of GW-BSE damped response calculations the targeted core states should be included in the orbital range when the GW quasiparticle energies are evaluated. A picture-change correction of the dipole moment is available for relativistic all-electron Hamiltonians ($rdkh, $rbss, $rx2c, also in their local variant $rlocal)and enforced by $pcc. For details on relativistic effects please, see Sec. 6.5.

8.4.4 Dynamic First Hyperpolarizability

Hyperpolarizability calculations are run through escf and can be requested by including

$scfinstab hyperpol <unit>
<freq1>
<freq2>
<pair1_a> <pair2_b>
...

in the control file. A hyperpolarizability calculation will be performed using the pairs given on each line AND using all combinations of the single frequencies given alone on each line. Specifying no frequencies initiates a static calculation. The same unit specifications are accepted as for dynamic polarizability calculations. For example,

$scfinstab hyperpol nm
800
1000
1064 1064

will compute hyperpolarizability tensors with frequency pairs

• 1000 nm, 1000 nm

• 1000 nm, 800 nm

• 800 nm, 800 nm

• 1064 nm, 1064 nm

8.4.5 Stability Analysis

Stability analysis of spin-restricted closed-shell ground states is enabled by

$scfinstab singlet

for singlet instabilities,

$scfinstab triplet

for triplet instabilities (most common), and

$scfinstab non-real

for non-real instabilities.

$scfinstab complex

for general complex instabilities (2c Kramers symmetric reference).

After that, it is necessary to specify the IRREPs of the electronic Hessian eigenvectors (“orbital rotations”) to be considered. Without additional knowledge of the system one usually needs to calculate the lowest eigenvalue within every IRREP:

$soes all 1

Positivity of the lowest eigenvalues in all IRREPs is sufficient for stability of the ground state solution. If one is interested in, say, the lowest eigenvalues in IRREPs e\(_{\text{g}}\) and t\(_{\text{2g}}\) only, one may specify:

$soes
eg   1
t2g  1

Triplet instabilities in the totally symmetric IRREP indicate open shell diradical states (singlet or triplet). In this case, start MOs for spin-symmetry broken UHF or UKS ground state calculation can be generated by specifying

$start vector generation

escf will provide the start MOs (\(\rightarrow\) $uhfmo_alpha, $uhfmo_beta) as well as occupation numbers (\(\rightarrow\) $alpha shells, $beta shells) for a spin-unrestricted calculation with equal numbers of \(\alpha\) and \(\beta\) electrons (pseudo-singlet occupation).

8.4.6 Vertical Excitation and CD Spectra

The calculation of excited states within the TDHF(RPA)/TDDFT approach is enabled by

$scfinstab rpas

for closed-shell singlet excitations,

$scfinstab rpat

for closed-shell triplet excitations, and

$scfinstab urpa

for excitations out of spin-unrestricted reference states.

If it is intended to use the TDA instead, specify

$scfinstab ciss

for closed-shell singlet excitations,

$scfinstab cist

for closed-shell triplet excitations,

$scfinstab ucis

for excitations out of spin-unrestricted reference states, and

$scfinstab spinflip

for spin-flip (\(z\)-component of the total spin changes by \(\pm 1\)) excitations out of spin-unrestricted reference states. For details concerning the theory see ref. [232]. In practice, this functionality can be used for the calculation of triplet-singlet, quartet-doublet, … excitations (see ref. [233] also for further information about the implementation). It is only available within the TDA in combination with LDA functionals and the HF exchange. It is strongly recommended to increase $escfiterlimit.

In the two-component case, specify

$scfinstab soghf

for two-component excitation energy calculations on closed-shell systems. [28] This implementation is only available in combination with LDA and GGA functionals; since version 7.4 also hybrid functionals are supported[30].

$scfinstab tdasoghf

for two-component excitation energy calculations on closed-shell systems using the TDA, where in addition HF exchange is accessible. [29]

open-shell systems can be accessed using TD-HF or the Bethe-Salpeter equation $bse

The keywords $soghf and $kramers in case of closed-shell systems also have to be set. Note that in two-component TD-DFT with metaGGAs and up, the current-dependent response is set as default. In two-component local hybrid functional TD-DFT calculations, the mixing between the exact exchange energy density and the LMF is neglected.

The Bethe-Salpeter equation (BSE) for excited states can be invoked by adding the keyword $bse. The correlation-augmented Bethe-Salpeter equation (cBSE) for excited states can be invoked by adding the keyword $cbse[31]. In BSE calculation RI-K is mandatory and automatically set.

Next, the IRREPs of the excitations need to be defined, which is again accomplished using $soes. For example, to calculate the 17 lowest excitations in IRREP b\(_{\text{1g}}\), the 23 lowest excitations in IRREP e\(_{\text{u}}\), and all excitations in IRREP t\(_{\text{2g}}\), use

$soes
b1g  17
eu   23
t2g  all

and run escf. Since point group symmetry cannot be exploited in two-component calculations, there is only the totally symmetric IRREP a.

Note that $soes specifies the IRREP of the excitation vector which is not necessarily identical to the IRREP of the excited state(s) involved. In general, the IRREP(s) of the excitation(s) from the ground to an excited state is given by the direct product of the IRREPs of the two states. For example, to calculate the first A\(_2\) state in a \(C_{\text{2v}}\)-symmetric molecule with a B\(_2\) (open-shell) ground state, it is necessary to specify

$soes
b1   1

The number of excitations that have to be calculated in order to cover a certain spectral range is often difficult to determine in advance. The total number of excitations within each IRREP as provided by the define ex menu may give some hint. A good strategy is to start with a smaller number of excitations and, if necessary, perform a second escf run on a larger number of states using the already converged excitation vectors as input.

To compute absorption and CD spectra, it is often sufficient to include optically allowed transitions only. This leads to substantial reduction of computational effort for molecules with higher symmetry. For example, in the UV-VIS spectrum of an O\(_{\text{h}}\) symmetric molecule, only \(t_{\text{1u}}\) excitations are optically allowed. The IRREPs of the electric and magnetic dipole moments as well as of the electric quadrupole moment are displayed automatically in the define ex menu.

If a large number of states is to be calculated, it is highly recommended to provide extra memory by specifying

$rpacor m

the integer m being the core memory size in megabytes (default is 20). The larger m, the more vectors can be processed simultaneously without re-calculation of integrals. As a rule of thumb, m should be ca. 90% of the available main memory. If RI-\(J\) is used ($ridft), it is recommended to set $ricore to a small value and $rpacor to a large value if the number of states is large, and vice versa if it is small. Since two-component calculations are more demanding concerning computation time and required memory it is strongly recommended to increase $rpacor.

By specifying

$spectrum unit

and/or

$cdspectrum unit

a list of excitation energies and oscillator and/or rotatory strengths of the optically allowed transitions is written onto file spectrum and/or cdspectrum. As above, unit specifies the energy unit and may be ev, nm, 1/cm, or a.u. (default). The files spectrum and cdspectrum may conveniently be used for further processing, e.g., using a plotting program such as Gnuplot.
Additionally, a spectrum broadened by gaussian functions can be generated by the Peak ANAlyzing MAchine (panama).[234] It reads in excitation energies and their corresponding oscillator strengths from the escf output file and prints the resulting spectrum to data.plot (on an eV axis) or to data.nm.plot (on a nanometer axis) which can be plotted by programs such as Gnuplot.
Both differential densities and non-relaxed difference densities can be visualized to get an impression of the character and localization of the excitation(s). Differential densities (ed.plt) are generated by egrad after setting the keyword $pointval in combination with the -proper option (see Sec. 22.2). A computational less demanding alternative are non-relaxed difference densities which are directly obtained from the escf output file by first running panama and subsequently dscf -proper or ridft -proper.[234]
Exact TDDFT transition density can also be plotted by the escf program, see section 22.2.

By specifying

$curswitchdisengage

inclusion of the current-density response for MGGA calculations is disabled. Note that the results of calculations using this flag will no longer be gauge-invariant and will differ from results obtained with the standard gauge-invariant implementation.
In general, CD spectra calculated with the length representation of the electric transition dipole moment are not gauge-invariant. Using GIAOs for the magnetic transition dipole moment, however, makes CD spectra gauge-invariant even if the length representation is used for the electric transition dipole moment. Calculating gauge-invariant rotatory strength tensors and thus CD spectra is possible for closed-shell molecules if the keyword $mgiao is added. First, an mpshift calculation needs to be run in order to obtain the perturbed density which is written on disk in a filed called ’umunu’. Then, a regular escf calculation can be run in order to calculate magnetic transition dipole moments which lead to gauge-invariant rotatory strength tensors for the length representation.

8.4.7 Two-photon absorption

2PA calculations are run through escf and can be requested with

$scfinstab twophoton <excited state method>
<freq1>
...
$soes
...
$exopt
...

where <excited state method> should be a description of the excited state method (rpas, ciss, urpa, ucis). A 2PA tensor will be computed for each excited state specified by the $soes block using the frequencies <freq1> and \(\Omega_n -\)<freq1>. Also, half may be provided as a frequency, which then uses half the excitation energy for each frequency (this is the default). 2PA amplitudes for specific states can be computed by specifying the states in $exopt. $exopt for different irreps can be specified just as for $soes or also as a comma separated list of indices (e.g., “3”) and ranges (e.g. “5-7”). For example

$scfinstab twophoton rpas
$soes
a1 6
a2 8
$exopt
a1 1-3, 5
a2 all

will compute

• 6 excited states in the a\(_1\) IRREP and 8 excited states in a\(_2\) and

• 2PA amplitudes with frequencies equal to \(\Omega_n/2\) for states 1, 2, 3, and 5 in irrep a\(_1\), and all 8 computed states in a\(_2\).

2PA amplitudes require the calculation of dynamic polarization vectors \(|X^{(\alpha)},Y^{(\alpha)}\rangle\) at the frequencies specified (e.g., \(\Omega_n/2\)). If 2PA amplitudes are requested for many states, these frequencies can become (nearly) resonant with lower-lying excitation energies, which can cause severe convergence problems. If such problems are encountered, try limiting the number of 2PA amplitudes computed in a single pass using $exopt.

8.4.8 Excited State Geometry Optimizations

The input for computing excited state gradients and properties using egrad is exactly the same as for an excited state calculation using escf, see the previous section. Gradients and properties are calculated only for one state at a time. By default, this is the highest excitation specified by $soes (only one IRREP is allowed). Sometimes, e.g. close to excited state intersections, it may be necessary to include higher excited states in the initial excitation vector calculation to prevent root flipping. This is accomplished using

$exopt n

which explicitly enforces treatment of the \(n\)-th state; n must be less or equal the number of states specified in $soes.

After the input for the ground and excited state calculations has been set up, an excited state geometry optimization can be started by issuing the command

nohup jobex -ex &

The option -ex forces jobex to call egrad instead of grad (or rdgrad if -ri is also specified). In each geometry step, the excitation energy is written on the fourth column in $energy, and the data group $last excitation energy change is updated. Otherwise, the excited state optimization proceeds in exactly the same way as a ground state optimization (see Chapter 3).

8.4.9 Excited State Force Constant Calculations

Excited state vibrational frequencies can be calculated by numerical differentiation of analytic gradients using NumForce (see Chapter 15). A NumForce calculation for an excited state may be started by the command

where n is the number of the excited state in \(C_1\) symmetry. In order to determine n, it is recommended to perform an escf calculation in \(C_1\) symmetry. Note that numerical calculation of excited state force constants is likely to fail if there are other states nearby (in \(C_1\)), because the roots may flip when the molecule is distorted. Note also that it may be necessary to include higher excited states (using $exopt, see above) in \(C_1\) calculations of molecules with higher symmetry in order to enforce convergence to the correct state. In any case, it should be checked that the energy change due to the displacements (available in the numforce/KraftWerk/*.log files) is reasonably small.

For a NumForce run, the convergence criteria should be tightened. It is recommended to use at least

$scfconv 8

in all NumForce calculations. Other NumForce options such as -central, -d, -np work in exactly the same way as they do for ground states.

8.4.10 Polarizability Derivatives and Raman Spectra

Calculations of polarizability derivatives by the egrad program use the same specifications in the $scfinstab data group as polarizability calculations by escf.

$scfinstab polly

specifies derivatives of the static polarizability, while

requests derivatives of the dynamical polarizability at the given frequency. Note that, unlike polarizability calculations, multiple frequencies are not allowed. Polarizability derivatives have to be projected onto vibrational normal modes to obtain Raman intensities, see Chapter 15 for further details.

8.4.11 State-to-state properties

All state-to-state properties, including transition moments and derivative couplings, are computed using the state-to-state derivative coupling machinery in egrad, as the state-to-state 1TDM is a byproduct of computing the derivative coupling. Only C\(_1\) symmetry is supported.

To trigger the calculation of state-to-state derivative couplings, an option must be given to the $nacme keyword.

$nacme <full/pseudo or response>

Providing either full or pseudo will compute derivative couplings under the pseudowavefunction approximation. This is the recommended option since pseudowavefunction couplings are well-behaved and stable. Providing response will compute derivative couplings with quadratic response theory. Couplings computed within response theory can diverge unphysically, so caution is advised.

For derivative couplings, it is also recommended to neglect the antisymmetric overlap integrals which are responsible for translational variance. This is approximately equivalent to incorporating electron translation factors (ETF) and is enabled by placing

$do_etf

in the control file.

The $coupled states keyword will control which states are coupled. $coupled states understands a comma separated list of numbers and ranges. Couplings will be computed between every pair of states specified on $coupled states. By default, $coupled states has the value “all”, which computes couplings between every state specified in $soes. For example,

$coupled states 1-3, 6, 8

will compute couplings between each pair of states in the list 1, 2, 3, 6, 8.

$coupled states all

will compute couplings between all available states.

When using $nacme is found, the $exopt key has added flexibility that allows the simultaneous calculation of multiple excited state gradients. It understands the same syntax as $coupled states such that

$coupled states 1-3, 6, 8

will compute gradients corresponding to states 1, 2, 3, 6, 8.

$coupled states all

will likewise compute gradients for all available states.

8.4.12 Nuclear spin-spin coupling constants

See Sec. 17.

8.4.13 Magnetic fields

TD-DFT and GW-BSE calculations in magnetic fields can be performed since V7.7. A two-component formalism is necessary for this, so $soghf needs to be set for magnetic fields. Note that for metaGGAs the full current-dependent kernel is used, effectively leading to current-dependent cTD-DFT.

8.4.14 Berry curvature

The Berry curvature tensor can be computed in a two-component ($soghf) framework: \[\begin{equation} \Omega_{I\alpha,J\beta} = 2 \text{Im} \ensuremath{\langle \frac{\partial \Psi}{\partial R_{I\alpha}} | \frac{\partial \Psi}{\partial R_{J\beta}} \rangle} \end{equation}\](8.27) This is an important quantity for the description of nuclear motion in the presence of external or internal magnetic fields. Furthermore, it is associated with a variety of non-adiabatic effects. In this implementation, the source of a non-vanishing Berry curvature can either be the presence of an external magnetic field ($magnetic field) or relativistic effects such as spin-orbit coupling (SOC, $rx2c). The Berry curvature can be calculated by using the following keywords:

$scfinstab polly
$berry

The $berry keywords prints out the Berry curvature at the end of an escf calculation, in addition to a variety of other related quantities. Furthermore, the Berry curvature tensor is printed into a file named ’berrycur’. Currently, the Berry curvature can be calculated on the Hartree-Fock (magnetic fields, SOC) and DFT (magnetic fields) level. The two-electron part has to be included without any approximations ($coulex).

The use of very tight convergence criteria (at least $denconv 1.0d-10 and $rpaconv 8) is recommended. Please note that the Berry curvature can be quite sensitive to the presence of instabilities in the wave function. If the calculation does not converge, check for instabilities ($scfinstab complex).

8.4.15 Diagonal Born-Oppenheimer correction

The diagonal Born-Oppenheimer correction (DBOC), \[\begin{equation} E^{DBOC} = \sum_I^{N_{\text{nuc}}} \frac{1}{2 M_I} [ \ensuremath{\langle \Psi | \hat{P}_I^2 \Psi \rangle} - \ensuremath{\langle \Psi | \hat{P}_I \Psi \rangle}^2] \, , \end{equation}\](8.28) is also calculated using the $berry keyword. It can be used as an estimate for non-adiabaticity in the framework of the adiabatic Born-Oppenheimer approximation. Since it is calculated alongside the Berry curvature, it is available using all methods described in the above section.

8.4.16 Prediciting colors using the Color Prediction Tool cpt

To obtain the calculated emission and absorption colors of a compound after a TD-DFT/GW-BSE calculation simply execute cpt in the same directory (basically only the exspectrum or spectrum file needs to be present). In case of response calculations from ricc2 or pnoccsd the $spectrum keywords needs to be present for the spectrum file to be generated. escf, starting from version 7.4, will always generate exspectrum files that can automatically be read in from cpt. Note: if more than one spectrum is present in the exspectrum/spectrum file cpt will generate the accumulated color. For some examples see test case TURBOTEST/escf/short/CPT and the included files and CRIT file.

8.4.17 Approximations for Coulomb and Exchange integrals

In hybrid functional TD-DFT (gradient) calculations most of the time is spent in evaluating the exchange-integrals. Therefore two ways to speed these up have been implemented:

$rick

use the RI-K approximation (recommended if enough RAM available)

$senex

use seminumerical exchange. The default settings of $esenex are recommended and lead to negligible errors in the excitation energies.

Usage of RI-K is highly recommended if enough RAM and a fast disk can be provided. RI-K is available for all types of calculations in the modules escf and egrad up to D2h symmetry. The keyword $rick is exclusive to escf, egrad and aoforce, other programs will not trigger usage of RI-K even if those parts will. The keyword $rik triggers RI-K in all modules where it is available, also (but not only) in escf, egrad and aoforce. RI-K is not available for range-separated hybrids. Seminumerical exchange is supported for all types of calculations in escf and egrad, also two-component calculations are fully supported. We recommend to use the default settings as they will yield reliable results in virtually all cases.

$esenex

Further, if one also wants to compute the Coulomb contribution using seminumerical techniques the $pseudospectral keyword may be added to the control file.

$esenex
$pseudospectral

Seminumerical exchange and fully pseudospectral approaches can be used with all point groups and functionals (global and range-separated hybrids). For a more detailed list of options please refer to the general keyword section.

8.4.18 Minimal auxiliary basis set

The calculation of TDDFT response properties can be speeded-up by about two-orders of magnitude using a minimal auxiliary basis set composed by just one s-type basis function per atom [235, 236] and switching off the calculation of the XC kernel on the grid. The resulting approach is named TDDFT-as[235], for semilocal XC functionals using the RI-J approach with a minimal auxiliary basis set for the Coulomb interaction, or TDDFT-ris[236], for hybrid XC functionals using the RI-K approach with the same minimal auxiliary basis set for both the Coulomb and exchange interaction. The resulting average deviation from conventional TDDFT excitation energies is less than 0.1 eV for organic molecules[236] and less than 0.02 eV metal clusters[235]. Note that we are considering only singlet excitation energies for closed-shell systems: open-shell systems and triplet excitations are inaccurate.

The TDDFT-as/TDDFT-ris approaches are based on the exact KS ground-states and approximate the linear-response. Thus these methods are similar to other tight-binding TDDFT methods[237, 238, 239] and less empirical. The exponents of the minimal auxiliary basis set need to be properly selected as they also approximate the XC kernel effects.

The keyword

$escfnoxc

disables the grid and the calculation of the XC kernel.

The script escfrisprep can be use to add this keyword to the control file as well as to include the minimal auxiliary basis set.

How to perform a TDDFT-as/TDDFT-ris:

• Run a conventional ground-state KS with RI-J or RI-JK approximation

• Use define to add all the required TDDFT keywords in the control file.

• Run escfrisprep

• Run escf

escfrisprep, with the default option (-m auto) will detect the XC type (reading the $xctype group) and define the appropriate basis set. Otherwise, for hybrid/RSH functionals, the TDDFT-ris method can be forced with:

escfrisprep -m ris

and the minimal basis set will be set in the $cbas group. For semilocal functionals, the TDDFT-as method can be forced with:

  
escfrisprep -m as 

and the minimal basis set will be set in the $jbas group.

For transition metal-complexes, the minimal auxiliary basis set is not sufficient for the description of low-lying excited states[236]. For those systems run:

escfrisprep -x element

where element is the element (e.g. the transition metal atom) to exclude. Use -x multiple times if you want to exclude more than one element (e.g -x Au -x Ag). For those atoms the full (default) auxiliary basis set will be used.

TDDFT-as/TDDFT-ris can also be used to generate initial guesses for further runs with standard TDDFT, which can reduce the number of iterations in some cases by a factor of  1.5. To do this, set up your calculation as previously described. Then run:

escfrisprep  
escf > escf-ris.out
escfrisprep -r 
escf > escf.out 

The command escfrisprep -r will restore the standard TDDFT files.