Table of Contents

Time-dependent density functional theory

Excitation energies and linear response properties can be calculated utilising the time-dependent density functional theory (TDDFT) method. The program should normally be called after a Kohn-Sham or Hartree-Fock calculation because it looks for the most recent orbital dump record to read in the MO coefficients and orbital energies. Further settings like the functional type and quadrature grid are then adopted from the previous ground-state calculation, yet some may also be modified via input commands in the TDDFT program. Detailed descriptions for the different calculation modes are given in the following subsections.

TDDFT program

Assuming that the molecule has $C_{2v}$ symmetry, a typical input for calculating the 6 lowest roots of the Hessian for IRREp 1 and 3, respectively, reads 

 ks,<func>
 tddft,states=[-6.1,-6.3]

The type of the xc kernel is usually adjusted automatically using the parameters from the previous ground-state KS calculation. However, various user inputs are available to change this, see below. Currently, density functionals of the (hybrid) LDA or GGA type as well as range-separated LDA/GGA functionals are supported in TDDFT. Note, however, that currently unsupported meta GGA’s (or special xc potentials for which the kernel can not (easily) be derived) may be combined with existing adiabatic LDA (ALDA) xc kernels.

Excitation energies for spin-unrestricted wave functions can be computed, too, using, e.g. 

 uhf
 tddft,states=[...]

for time-dependent Hartree-Fock (TDHF) or 

 uks,<func>
 tddft,states=[...]

for time-dependent DFT.

The program can be run in various integral transformation modes using 

  tddft,inttype=<ityp>

with ityp=[ao|mo|mo(df)|df]. With ityp=ao (set as the default) an MO-AO transformation of the excitation vectors is performed. This requires the 4-indexed Coulomb repulsion integrals in the AO basis to be stored on disk. ityp=mo means that a full 4-index MO transformation to the Coulomb and exchange operators (2-external MO integrals) is performed. These are then used to build the Hessian-vector products in each TDDFT Davidson iteration step. While this option is faster than ityp=ao it requires a larger amount of disk space available.

While both, ityp=ao and ityp=mo require the storage of the 4-indexed AO integrals on disk, it is also possible to avoid this using density-fitting 

  df-tddft,...

in which case only 3-indexed integrals are used to compute the Hessian-vector products. While this requires an additional auxiliary basis set to be defined, the program can usually detect the fitting basis set that corresponds to the given orbital basis automatically, so that normally no additional information needs to be passed to the program. If a special fitting basis set shall be used one can add 

  df-tddft,basis=<auxbas>...

with auxbas corresponding to a basis set keyword that was defined in the basis section.

Note that density-fitting can also be used in conjunction with inttype=mo in which case the Coulomb and exchange integrals are computed using density-fitting, thus avoiding the storage of the 4-indexed AO integrals.

The following list summarises the keywords available in the TDDFT program:

Parameters which influence the behaviour of the Davidson solvers that can be used (should normally only be modified if convergence is hampered for some reason):

A number of options from the list above can be given separately in the TDDFT command group, e.g.: 

 {tddft,...; core,...; fxks,...}

Some of them are exclusive, namely:

Default values for all parameters may also be looked up in the tddft.registry file.

Analysing the spectrum

By default, Molpro prints oscillator strengths and the transition moments correponding to the excitations computed at the end of the calculation in a summary table. Using 

tddft,...; gnuplot,<spec.gp>

where spec.gp is a generic file name, the spectrum can be exported to a file which can be loaded with Gnuplot. For example, for the water molecule (PBE0 functional) the resulting plot for the absorption spectrum is h2o-spectrum.pdf. Note that the shape of the plot can be adapted by modifying the value of sigma in 'spec.gp which denotes the standard deviation in units of 1 nm$^{-1}$ for the Gaussian approximations for the respective transitions.

The excitation vectors can be exported to formats for visualisation by using either 

tddft,...; molden,<vecs.mold>

which exports the excitation densities to a file which can be read with the program Molden (vectors are stored as MO coefficients, so can be visualised using the orbital selection interface of Molden) or, alternatively, via 

tddft,...; cube,<vecs.cube>

in which case for each excitation a Gaussian cube file is created which can be visualised with various external programs, e.g., JMol. Note that these two options currently are implemented only for closed-shell restricted calculations.

For 'large' molecules one should be aware of the fact the TDDFT excitations can be strongly shifted compared to the experiment due to the charge transfer problem (see, e.g., https://doi.org/10.1103/PhysRevA.80.012507). It is strongly recommended to use hybrid or (much better) range-separated density functionals to obtain reliable results for the transition energies for such systems. Note that the failure for Rydberg excitations for standard DFT methods can be resolved using asymptotically corrected xc potentials, see section Asymptotic correction for xc-potentials(ASYMP) and the paper by Handy and Tozer http://dx.doi.org/10.1063/1.477711 for more details.

In order to be able to distinguish between the different types of excitations, apart from the plotting options described above, one can print the values of $\Lambda$ described in the paper by Peach, Benfield, Helgaker and Tozer https://doi.org/10.1063/1.2831900: 

 tddft,diagnost=1,...

The values are defined as $$\Lambda=\frac{\sum_{ia}\kappa_{ia}O_{ia}}{\sum_{ia} \kappa_{ia}}$$ with $$O_{ia}=\langle|\phi_i||\phi_a|\rangle=\int d{\bf r} |\phi_i({\bf r})||\phi_a({\bf r})|$$ where $\kappa_{ia}$ is a coefficient of the excitation vector and $O_{ia}$ corresponds to the spatial overlap between an occupied orbital $\phi_i$ and an unoccupied orbital $\phi_a$, see https://doi.org/10.1063/1.2831900 for more details. Note that the value of $\Lambda$ can depend on both the functional type as well as the basis set. As a rule of thumb, values of $\Lambda<0.3$ obtained with the B3LYP functional indicate a strongly nonlocal character (the smaller the more nonlocal)

Reponse properties from the coupled perturbed Kohn-Sham method

The linear response to (frequency-dependent) perturbations can be calculated with the TDDFT program using (use states=[] explicitly in order to skip the calculation of excitation energies): 

 tddft,states=[],...; prop,<op1>,<op2>,...; om,<om1>,<om2>,...

where <op> can be any operator given in the section One-electron operators and expectation values (GEXPEC). A list of frequencies ($\omega$) can be specified by om,… where positive values denote real frequencies while in case of negative values the response to perturbations oscillating at imaginary frequencies is calculated. Real and imaginary input values can be mixed arbitrarily. If static response properties are requested, too, the value of om=0.0 should be inserted at the beginning of the list (note that if the list of $\omega$'s is omitted in the input, the static response is calculated by default).

For the calculation of multipole-multipole polarisabilities it is possible to use the short-hand input variant 

 tddft,states=[],...; prop,q1,q2,q3

in which case all $l=1,2,3$ rank responses (dipole, quadrupole and octopole) are computed without having to insert all individual cartesian components. If all components for a given rank are given, the program performs a transformation to the corresponding spherical harmonics representation and prints the results in the output as well.

Dispersion coefficients

Dispersion coefficients can be calculated using 

  tddft,states=[],...; prop,...; disp,<nfreq>

where nfreq denotes the number of quadrature points for the numerical calculation of the integral over imaginary frequencies with the Casimir-Polder integral transform. Normally, values of the order of nfreq=$10$ are sufficient for accurate integrations. The Gauss-Legendre quadrature scheme as described in the paper by Amos et al. is used in the response program.

The isotropic leading order $C_6$ dispersion coefficient between two monomers $A$ and $B$ is given by $$C_6^{AB} = \frac{3}{\pi}\int_0^\infty d{\omega}~\alpha^A(i\omega)\alpha^B(i\omega) $$ with $\alpha^A$ and $\alpha^B$ denoting the isotropic dipole-dipole polarisabilities of the two monomers. The $C_6$ coefficients are computed automatically if prop,q1 is used. Higher order coefficients $C_8$ and $C_{10}$ can be computed as well if quadrupole-quadrupole and dipole-octopole polarisabilities are available (prop,q1,q2 or prop,q1,q2,q3). Note that the $C_{10}$ coefficients are computed even if no octopole moments are given in the input file. The values then contain only the quadrupole-quadrupole $\times$ quadrupole-quadrupole polarisability contributions.

For obtaining accurate results: use asymptotically corrected xc potentials 'and' basis sets with 'enough' diffuse functions!

The following options can be set in the TDDFT response program:

Rotatory strengths

Rotatory strengths for the excitations requested can be calculated with 

{tddft,states=...,cdspec=1,...}

These can be used to simulate CD spectra for optically active molecules. Since the conventional length representation of the underlying transition moments is not gauge-invariant (and so can strongly deviate from the cbs result for small basis sets), also the velocity representation for the electric transition dipole moments are computed and printed in the output, see, e.g., works by Autschbach et al. https://doi.org/10.1063/1.1436466.

Rotatory strenghths can be currently computed only for the spin-restricted TDDFT case.

Excited state gradient calculations

For performing TDDFT gradient calculations and geometry optimisations for excited states the value of the grad parameter needs to be set to 1 (>0) in the tddft command: 

    tddft,grad=1,....

Currently this is possible for closed-shell singlet excited states using LDA, GGA, hybrid-GGA or range-separated hybrid functionals.

In TDDFT gradient calculations it is recommended to set tighter thresholds for the convergence of the excited state vectors and for the iterative solution of the Z-vector equation. E.g., for small systems the following settings would be reasonable: 

    tddft,grad=1,thrv=1d-8,thrcg=1d-12,....

while for larger systems the threshold values may be set larger but should normally be more tight compared to conventional TDDFT energy calculations, see above. The TDDFT program then computes the excited state density matrix as well as the excited state contribution to the overlap lagrangian and writes these to internal records for later use. The gradient can then be computed with the force program, see chapter Energy gradients for more details: 

    tddft,grad=1,....
    force

and geometry optimisations can be done with the optg command: 

    tddft,grad=1,....
    optg

Note that by default the gradient is computed for the lowest excitation of a given symmetry (calculations for multiple symmetries are not allowed for gradient calculations). This can be changed by setting the parameter gradvec, e.g. 

    tddft,grad=1,states=[-4.2],gradvec=4,....

will compute the gradient for the 4th excited state of the second symmetry. It is necessary, of course, that the number of states requested by the states command is equal to or larger than gradvec. Note that an eigenvector-following method like is implemented in the EOM-CCSD program is not yet available for TDDFT geometry optimisations. Because of this it might happen that the initial starting guess is not followed strictly when there are crossings with other close excitations on the PES. This will be improved in a future version of the program.

There are a few special parameters which can influence the behaviour of tddft gradient calculations. Most of them are for debugging purpose and should normally not changed compared to their default values set in the registry file.

Nonadiabatic coupling matrix elements

First order nonadiabatic coupling matrix elements (NACME's) between the ground state and an excited state can be computed using TDDFT response theory (and so can be obtained even when the wave functions of the respective states are not known explicitly). In general a NACME between two states $n$ and $m$ is given as $$\tau^\xi_{nm}=\langle \Psi_n|\frac{\partial}{\partial \xi}|\Psi_m\rangle $$ where $\xi$ is a cartesian nuclear coordinate and $\Psi_n$ and $\Psi_m$ are the wave functions of the two states. In order to derive an expression for the first-order NACME from time-dependent reponse theory one considers the time evolution of the imaginary matrix element $$C_\lambda^\xi(t)=\langle \Psi_\lambda(t)|\frac{\partial}{\partial\xi}|\Psi_\lambda(t)\rangle$$ under the influence of a monochromatic one-particle perturbation $\hat{W}$ and finds that the NACME is given by the residues of $\overline{C}^{\xi(1)}$ at the excitation energies $\omega_n$ $$\tau_{0n}^\xi = \frac{2\pi i\mbox{Res}[\overline{C}^{\xi(1)}(\omega); \omega_n] }{\langle \Psi_0|\hat{\overline{W}}(\omega)|\Psi_n\rangle}$$

Following the work of Send and Furche https://doi.org/10.1063/1.3292571 this term can be computed using TDDFT ground-state response theory, thus without knowledge of the excited state wave function. The resulting expression that was derived by Send and Furche includes Pulay type terms that greatly reduce the basis set dependence of the NACME's and that reduces to the exact Chernyak–Mukamel formula https://doi.org/10.1063/1.480511 for first-order NACMEs in the complete basis-set limit, see https://doi.org/10.1063/1.3292571.

In Molpro, first-order NACME's can be computed with the TDDFT program by setting the logical-type option nacme=1 (true), for example: 

{tddft,states=[-4.1],nacme=1,nacmevec=2,...}

The option nacmevec can be used to set the state $n$ for which $\tau_{0n}^\xi$ shall be computed (otherwise NACME's for all states requested will be calculated). In the above example the four lowest excitations are computed and then the NACME $\tau_{02}^\xi$ is computed. Of course, nacmevec must be lower than or equal to the number of excitations requested.

First-order NACME's from TDDFT response theory have been implemented for closed and open-shell Hartree-Fock, LDA and GGA-type functionals.