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).
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,
or
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.
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 400 i nm (i is the imaginary unit), specify
.
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}ω^{2}∕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
Stability analysis of spin-restricted closed-shell ground states is enabled by
for singlet instabilities,
for triplet instabilities (most common), and
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).
The calculation of excited states within the TDHF(RPA)/TDDFT approach is enabled by
for closed-shell singlet excitations,
for closed-shell triplet excitations, and
for excitations out of spin-unrestricted reference states.
If it is intended to use the TDA instead, specify
for closed-shell singlet excitations,
for closed-shell triplet excitations,
for excitations out of spin-unrestricted reference states, and
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
for two-component excitation energy calculations on closed-shell systems. [121] This
implementation is only available in combination with LDA and GGA functionals.
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_{2} state in a C_{2v}-symmetric molecule with a B_{2} (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
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
and/or
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
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.
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
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).
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
nohup NumForce -ex n > force.out &
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
in all Numforce calculations. Other Numforce options such as -central, -d, -np work in exactly the same way as they do for ground states.
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
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.