## 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:

A list which supplies the number of excitations to be calculated per IRREp: [N$_1$.IRREP$_1$,N$_2$.IRREP$_2$,$\dots$]`STATES`

Sets the auxiliary basis set if`BASIS`

`MODE`

=1,`DFIT=.true.`

Factor that determines the amount of the ALDA exchange (Slater-Dirac) kernel (if`FXKS`

`XCMODE`

$\le 10$).Factor that determines the amount of the ALDA correlation (VWN5) kernel (if`FCKS`

`XCMODE`

$\le 10$).Factor that determines the amount of exact Hartree-Fock exchange.`FXHF`

XC kernel type whuch can take the settings:`XCKERNEL`

`XCKERNEL=[lda|gga|lda(df)|gga(df)]`

. Concerning`lda(df)|gga(df)`

, see`FITFXC`

option below.Sets the type of exchange-correlation kernel routines to be used in the Hessian-vector transformations.`XCMODE`

`XCMODE`

=1-10 $\rightarrow$ reserved for ALDA kernels,`XCMODE`

=11-20 $\rightarrow$ reserved for GGA kernels. (normally not required)Set to $\ne 0$ if triplet excitations instead of singlet excitations shall be computed.`TRIPLET`

Set the number of core orbitals.`CORE`

Threshold for orbital screening (default: $1e-5$).`TOLORB`

Block size of numerical quadrature batches used in xc kernel integrations.`NBLOCK`

Compute the xc kernel contribution to the Hessian via density-fitting in df-tddft calculations. With`FITFXC`

`FITFXC`

=$1$ the xc kernel integrals in the auxiliary basis are computed in advance and then are simply loaded during the DF-construction of the Hessian-vector products, which is usually the fastest option. With`FITFXC`

=$2$ the xc contribution to the Hessian is fitted via a contraction with the density fitting coeffiencients during the Davidson iterations. (default:`FITFXC`

=$0$)Use a singularity correction to fix the potentially singular auxiliary integrals over the xc kernel, see https://doi.org/10.1063/1.4893990. (default:`SING`

`SING`

=1)Parameter for singularity correction in FXCFIT case (see https://doi.org/10.1063/1.4893990. default:$1d-5$).`ESHIFT`

Parameter for singularity correction in FXCFIT case (see https://doi.org/10.1063/1.4893990. default:$1d-4$).`FSCAL`

Skip quadrature points around dummy centres if the density is lower than the given threshold when`THRRHODUM`

`FITFXC`

$>0$ is used. (default:`THRRHODUM`

=0d0)Switch between dfit (0) and ri technique (1) for`FIT1MOD`

`FITFXC`

$=1$ case. (default:`FIT1MOD`

=$1$).If $>0$ compute the Coulomb kernel contribution via numerical quadrature rather than calculating the Coulomb integrals explicitly. Note: this option is not available with all settings and is normally only used in conjunction with range-separated density functionals in which case the full-range Coulomb integrals are not available on disk (`COULGRID`

`INTTYPE`

=$0,1$ case). (default:`COULGRID`

=$0$)Compute triplet excitation energies instead of singlet excitations (note: singlet and triplet excitations are obtained both simultaneously in case of unrestricted calculations, so the option only makes sense for the closed-shell case) (default:`TRIPLET`

`TRIPLET`

=$0$).Calculate the 4-indexed Fxc operator integrals if`FXCOP`

`inttype=mo(df)`

(default:`FXCOP`

=$0$)Perform a full diagonalisation of the Hessian if`FULLDIAG`

`inttype=mo(df)`

(default:`FULLDIAG`

=$0$)Print the $\Lambda$ values from Tozers diagnostic test in the summary table, see https://doi.org/10.1063/1.2831900 for details. This can help to identify local vs. Rydberg vs. nonlocal charge-transfer excitations. (default:`DIAGNOST`

`DIAGNOST`

=$0$)If symmetry is used, print a summary table including the excitation energies for each IRREp. (default:`PRINTSTEP`

`PRINTSTEP`

=0)Number of excitation vector coefficients to be printed in the summary table (default:`NCOEFFPRI`

`NCOEFFPRI`

=$4$)

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):

Maximum number of iterations (default: $50$).`MAXITDAV:MAXIT`

Threshold to test degenerate states. If, via the number of states requested, one 'cuts' through a degeneracy (as is tested through the orbital energy differences) the program adds further states to be computed in order to improve/enable convergency (default: $1d-5$).`DEGTHR:THRDEG`

Davidson solver mode that can have values 1-5. Default:`DAVMOD`

`DAVMOD`

=5 (recommended).Maximum number of expansion vectors in the eigensolver procedure.`MAXVEC`

Number of start unit vectors (default:`STARTVEC`

`STARTVEC`

=“0”, meaning that the program usually sets the number of vectors to`NSTATES`

at the beginning).Energy threshold value. Default:`THRE`

`THRE`

=1d-6.Residual threshold value (Euklidean norm of V). Default:`THRV`

`THRV`

=1d-4.Threshold for smallest eigenvalue of S (`THRS`

`DAVMOD`

=5, default:`THRS`

=1d-6).Threshold for smallest eigenvalue of S of update vectors (`THRS_UPDATE`

`DAVMOD`

=5, default:`THRS`

=1d-4).Threshold for norm of expansion vectors (`THRD`

`DAVMOD`

=5, default:`THRS`

=-1d0).If the norm of a new computed expansion vector is smaller than this value, it will be discarded in the iterative update of the vector space (for`THRQ`

`DAVMOD`

=1,2). Default:`THRR`

=1d-2.If`USYM`

`DAVMOD`

=1,2 either compute a symmetrised subspace Hessian (`USYM`

=0) or not (`USYM`

=1). Default:`USYM`

=0.Perform deflation step in the Davidson solver (`DEFL`

`DAVMOD`

=1,2). Default:`DEFL`

=1.Denominator shift (`SHIFT`

`DAVMOD`

=5, default:`SHIFT`

=0d0).If true, project new vectors against old ones`PROJECT`

`DAVMOD`

=5, default:`PROJECT`

=1).If $>0$ print energy checks at the end (`CHECK`

`DAVMOD`

=5, default:`CHECK`

=0).If true, use exact hessian diag for lowest vectors (`REPLACE_DIAG`

`DAVMOD`

=5, default:`REPLACE_DIAG`

=0).If 0, compute hred directly; if 1, hred=bred*ared (`SOLVER_MODE`

`DAVMOD`

=5, default:`SOLVER_MODE`

=-1).

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:

Perform a stability analysis of the SCF wave function (computes the lowest eigenvalue of the electronic Hessian ${\bf A}+{\bf B}$).`STABIL`

Save the excitation vectors to a record for, e.g., restarting another TDDFT calculation from these`SAVE`

Read excitation vectors from a record and use these as the initial start vectors (it goes without saying that the dimensions etc. should match exactly those of the current calculation)`READ`

Truncate the number of virtual orbitals (only possible with`NVIR`

`inttype=mo(df)`

)Specify the number of core orbitals set exclusively for the exchange contribution to the Hessian, ie., this setting enables to use different sizes for the core for the Coulomb and the exchange contributions.`NCOREX`

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:

In case of unsymmetric Hessians (hybrid or range-separated functionals) one can choose between two different linear solvers (default:`LESOLVE`

`LESOLVE`

=$1$). Otherwise, a standard conjugate gradient solver is used to solve the response equations.Convergence threshold value for the (conjugate gradient) solver (default:`THRCG`

`THRCG`

=$1d-8$)

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

If set to`GRADEXC`

`GRADEXC`

=0 the ground-state DFT gradient contribution for the xc energy are computed with the standard DFT routines. When set to`GRADEXC`

=1 both the xc energy and kernel gradient contributions are computed simultaneously in a separate routine.Print the gradient at various stages in the program if set >0.`GRADPRI`

Use singular value decomposed density matrices in some AO based xc routines.`DENSVD`

Switch between different df-gradient modes for the calculation of the 2el-contribution to the tddft gradient.`DFGRADMODE`

Treshold for singular value decompositions of AO-basis density matrices.`THRDDEC`

Maximum number of auxiliary functions for transformed integral blocks in the density-fitting Fock routine.`MAXAUX`

Maximum number of loop batches in df gradient routines.`MAXBAT`

Use disk (=1) or a global array (=0) to store half-transformed integrals in the density-fitting Fock routine. This setting will only be used during the calculation of the gradient terms.`DISK`

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