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

in control, perform a dscf statistics run, if semi-direct integral processing is to be used (see Chapter 3.1), 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 increase the default grid size for the evaluation of the DFT functional from m3 to m4. In rare cases avoiding multiple grids can also increase numerical stability (use 4 instead of m4). To perform a two-component TDDFT calculation, the two-component version of ridft has to be run before (see Chapter 6.4) using the keywords $soghf and $kramers.

8.4.2 Polarizabilities and Optical Rotations

The calculation of dynamic polarizabilities is controlled by the keyword

$scfinstab dynpol unit
list of frequencies

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 400 i nm (i is the imaginary unit), specify

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 [α]_ω in deg⋅[dm⋅(g/cc)]^-1 is given in Eq. (15) of ref.  [119].

(8.17)

where C = 1.343 ⋅ 10^-4ω²∕M with M being the the molar mass in g/mol, ω the frequency in cm^-1, and δ(ω) is 1/3 trace of the electronic rotatory dispersion tensor given in atomic units.
Please note, that δ(ω) 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

before starting escf.

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

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:

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 eg and t2g only, one may specify:

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

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

8.4.4 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 ±1) excitations out of spin-unrestricted reference states. For details concerning the theory see ref. [126]. In practice, this functionality can be used for the calculation of triplet-singlet, quartet-doublet, ... excitations (see ref. [127] 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. [121] This implementation is only available in combination with LDA and GGA functionals.

$scfinstab tdasoghf

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

The keywords $soghf and $kramers also have to be 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 b1g, the 23 lowest excitations in IRREP eu, and all excitations in IRREP t2g, use

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₂ state in a C_2v-symmetric molecule with a B₂ (open-shell) ground state, it is necessary to specify

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_h symmetric molecule, only t_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). [128] It reads in excitation energies (in eV) and their corresponding oscillator strengths from the escf output file and prints the resulting spectrum to data.plot 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. 18.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. [128]

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.

8.4.5 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

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

8.4.6 Excited State Force Constant Calculations

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

where n is the number of the excited state in C₁ symmetry. In order to determine n, it is recommended to perform an escf calculation in C₁ symmetry. Note that numerical calculation of excited state force constants is likely to fail if there are other states nearby (in C₁), 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₁ 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

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

specifies derivatives of the static polarizability, while

$scfinstab dynpol unit
frequency

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 14 for further details.