This is a short tutorial on how to run linear-response TDDFT computations for absorption and emission spectroscopy. The TDDFT module enables a description of excitation energies and excited-state computations within the Tamm-Dancoff approximation (TDA) featuring GGA and hybrid functionals as well as semi-empirical simplified TDA kernels. The details of the implementation can be found in https://aip.scitation.org/doi/full/10.1063/1.5078682 (Ref. [1]) and in https://pubs.acs.org/doi/10.1021/acs.jctc.2c00144 for corresponding excited-state gradients. Note that the current module is based on an earlier TDDFT implementation https://chimia.ch/chimia/article/view/2005_499. Please cite these papers if you were to use the TDDFT module for the computation of excitation energies or excited-state gradients.

The implementation in CP2K is based on the Tamm-Dancoff approximation (TDA), which describes each excited state $p$ with the excitation energy $\Omega^p$ and the corresponding excited-state eigenvectors $\mathbf{X}^p$ as an Hermitian eigenvalue problem

$$\begin{equation} \label{tda_equation} \begin{aligned} \mathbf{A} \mathbf{X}^p = \Omega^p \mathbf{X}^p \, , \\ \sum_{\kappa k} [ F_{\mu \kappa \sigma} \delta_{ik} - F_{ik \sigma} S_{\mu \kappa} ] X^p_{\kappa k \sigma} + \sum_{\lambda} K_{\mu \lambda \sigma} [\mathbf{X}^p] C_{\lambda i \sigma} = \sum_{\kappa} \Omega^p S_{\mu \kappa} X^p_{\kappa i \sigma} \, . \end{aligned} \end{equation}$$

The Hermitian matrix $\mathbf{A}$ contains as zeroth-order contributions the difference in the Kohn-Sham (KS) orbital energies $\mathbf{F}$, and to first order kernel contributions $\mathbf{K}$ which comprise –depending on the chosen density functional approximation – Coulomb $\mathbf{J}$ and exact exchange $\mathbf{K}^{\rm{\tiny{EX}}}$ contributions as well as contributions due to the exchange-correlation (XC) potential $\mathbf{V}^{\rm{\tiny{XC}}}$ and kernel $\mathbf{f}^{\rm{\tiny{XC}}}$,

$$\begin{equation} \label{f_and_k_matrix} \begin{aligned} F_{\mu \nu \sigma} [\mathbf{D}] &= h_{\mu \nu} + J_{\mu \nu \sigma} [\mathbf{D}] - a_{\rm{\tiny{EX}}}K^{\rm{\tiny{EX}}}_{\mu \nu \sigma} [\mathbf{D}] + V_{\mu \nu \sigma}^{\rm{\tiny{XC}}} \, , \\ K_{\mu \nu \sigma} [\mathbf{D}^{\rm{\tiny{X}}}] &= J_{\mu \nu \sigma} [\mathbf{D}^{\rm{\tiny{X}}}] - a_{\rm{\tiny{EX}}} K^{\rm{\tiny{EX}}}_{\mu \nu \sigma}[\mathbf{D}^{\rm{\tiny{X}}}] + \sum_{\kappa \lambda \sigma'} f^{\rm{\tiny{XC}}}_{\mu \nu \sigma,\kappa \lambda \sigma'} D_{\kappa \lambda \sigma'}^{\rm{\tiny{X}}} \, . \end{aligned} \end{equation}$$

$\mathbf{S}$ denotes the atomic-orbital overlap matrix, $\mathbf{C}$ the occupied ground-state KS orbitals and $\mathbf{D}$ and $\mathbf{D}^{\rm{\tiny{X}}}$ ground-state and response density matrices,

$$\begin{equation}\label{density_matrices} \begin{aligned} D_{\mu \nu \sigma} = \sum_k C_{\mu k \sigma} C_{\nu k \sigma}^{\rm{T}} \, , \\ D_{\mu \nu \sigma}^{\rm{\tiny{X}}} = \frac{1}{2} \sum_{k} ( X^p_{\mu k \sigma} C_{\nu k \sigma}^{\rm{T}} + C_{\mu k \sigma} (X^p_{\nu k \sigma})^{\rm{T}} ) \, . \end{aligned} \end{equation}$$

Within the current implementation, symmetrization and orthogonalization of the response density matrix is ensured at each step of the Davidson algorithm. The current implementation features to approximate the exact exchange contribution of hybrid functionals using the auxiliary density matrix method (ADMM). Furthermore, the standard kernel can be approximated using the semi-empirical simplified Tamm-Dancoff approximation (sTDA), neglecting in this case XC contributions and approximating both Coulomb and exchange contributions $\mathbf{J}$ and $\mathbf{K}$ using semi-empirical operators $\boldsymbol{\gamma}^{\rm{\tiny{J}}}$ and $\boldsymbol{\gamma}^{\rm{\tiny{K}}}$ depending on the interatomic distance $R_{AB}$ of atoms $A$ and $B$,

$$\begin{equation} \label{stda_kernel} \begin{aligned} \gamma^{\rm{\tiny{J}}}(A,B) &= \left ( \frac{1}{(R_{AB})^{\alpha} + \eta^{-\alpha}} \right)^{1/\alpha} \, , \\ \gamma^{\rm{\tiny{K}}}(A,B) &= \left ( \frac{1}{(R_{AB})^{\beta} +( a_{\rm{\tiny{EX}}}\eta)^{- \beta} } \right )^{1/\beta} \, , \end{aligned} \end{equation}$$

that depend on the chemical hardness $\eta$, the Fock-exchange mixing parameter $a_{\rm{\tiny{EX}}}$ and powers of $\alpha$ and $\beta$ for either Coulomb and exchange interactions.

Within the current implementation, oscillator strengths can be calculated for molecular systems in the length form and for periodic systems using the velocity form (see Ref.[1]).

Based on Eq.\ ($\ref{tda_equation}$), excited-state gradients can be formulated based on a variational Lagrangian for each excited state $p$,

$$\begin{equation} \begin{aligned} L [\mathbf{X}, \mathbf{C}, \Omega, \bar{\mathbf{W}}^{\rm{\tiny{X}}}, \bar{\mathbf{Z}}, \bar{\mathbf{W}}^{\rm{\tiny{C}}} ] &= \Omega - \sum_{\kappa \lambda k l \sigma} \Omega ( X_{\kappa k \sigma }^{\rm{T}} S_{\kappa \lambda } X_{\lambda l \sigma} - \delta_{kl} ) \\ &- \sum_{kl \sigma} ( \bar{W}_{kl \sigma}^{\rm{\tiny{X}}} )^{\rm{T}} \sum_{\kappa \lambda} \frac{1}{2} ( C_{\kappa k \sigma}^{\rm{T}} S_{\kappa \lambda} X_{\lambda l \sigma} + X_{\kappa k \sigma}^{\rm{T}}S_{\kappa \lambda} C_{\lambda l \sigma}) \\ &+ \sum_{\kappa k \sigma}( \bar{Z}_{\kappa k \sigma})^{\rm{T}} \sum_{\lambda} ( F_{\kappa \lambda \sigma}C_{\lambda k \sigma} - S_{\kappa \lambda } C_{\lambda k \sigma} \varepsilon_{k \sigma}) \\ &- \sum_{kl\sigma} (\bar{W}^{\rm{\tiny{C}}}_{kl \sigma})^{\rm{T}} ( S_{kl \sigma} - \delta_{kl})\, . \end{aligned} \end{equation}$$

introducing Lagrange multipliers $\bar{\mathbf{W}}^{\rm{\tiny{X}}}$, $\bar{\mathbf{W}}^{\rm{\tiny{C}}}$, and $\bar{\mathbf{Z}}$ to ensure stationarity of the corresponding ground-state (GS) equations and to account for the geometric dependence of the Gaussian orbitals and thus requiring to solve the Z vector equation iteratively.

To compute absorption spectra, parameters defining the LR-TDDFT computation have to be specified in the TDDFPT subsection of the section PROPERTIES of section FORCE_EVAL. Furthermore, RUN_TYPE has to be set to ENERGY and the underlying KS ground-state reference has to be specified in the DFT section.

The most important keywords and subsections of TDDFPT are:

• KERNEL: option for the kernel matrix $\mathbf{K}$ to choose between the full kernel for GGA or hybrid functionals and the simplified TDA kernel
• NSTATES: number of excitation energies to be computed
• CONVERGENCE: threshold for the convergence of the Davidson algorithm
• RKS_TRIPLETS: option to switch from the default computation of singlet excitation energies to triplet excitation energies

To compute excited-state gradients and thus corresponding fluorescence spectra, the excited state to be optimized furthermore has to be specified by adding the subsection EXCITED_STATES of the section DFT.

acetone_tddft.inp

&GLOBAL
  PROJECT S20Acetone
  RUN_TYPE ENERGY
  PREFERRED_DIAG_LIBRARY SL
  PRINT_LEVEL medium
&END GLOBAL
&FORCE_EVAL
  METHOD Quickstep
&PROPERTIES
  &TDDFPT                          ! input section for TDDFPT
   KERNEL FULL                       ! specification of the underlying kernel matrix K
                                     ! FULL kernel is for GGA and hybrid functional computations
                                     ! sTDA kernel is referring to a semi-empirical sTDA computation
   NSTATES 10                      ! specifies the number of excited states to be computed
   MAX_ITER   100                  ! number of iterations for the Davidson algorithm
   CONVERGENCE [eV] 1.0e-7         ! convergence threshold in eV
   RKS_TRIPLETS F                  ! Keyword to choose between singlet and triplet excitations
 !  &XC                            ! If choosing kernel FULL, the underlying functional can be
 !   &XC_FUNCTIONAL PBE0             ! specified by adding an XC section
 !   &END XC_FUNCTIONAL              ! If not choosing ADMM, the functional can be chosen
 !  &END XC                        ! independently from the chosen GS functional
  &END TDDFPT
&END PROPERTIES
  &DFT
    &QS
      METHOD GPW
      EPS_DEFAULT 1.0E-17
      EPS_PGF_ORB 1.0E-20
    &END QS
    &SCF
      SCF_GUESS restart
      &OT
         PRECONDITIONER FULL_ALL
         MINIMIZER DIIS
      &END OT
      &OUTER_SCF
         MAX_SCF 900
         EPS_SCF 1.0E-7
      &END OUTER_SCF
      MAX_SCF 10
      EPS_SCF 1.0E-7
    &END SCF
    POTENTIAL_FILE_NAME POTENTIAL_UZH
    BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
    BASIS_SET_FILE_NAME BASIS_ADMM_UZH
    &MGRID
      CUTOFF 800
      REL_CUTOFF 80
    &END MGRID
    &AUXILIARY_DENSITY_MATRIX_METHOD
      METHOD BASIS_PROJECTION
      EXCH_SCALING_MODEL NONE
      EXCH_CORRECTION_FUNC NONE
      ADMM_PURIFICATION_METHOD NONE
    &END AUXILIARY_DENSITY_MATRIX_METHOD
    &POISSON
       PERIODIC NONE
       POISSON_SOLVER WAVELET
    &END
    &XC
     &XC_FUNCTIONAL PBE0
     &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC [angstrom] 14.0 14.0 14.0
      PERIODIC NONE
    &END CELL
    &COORD
       C 0.000000 1.282877 -0.611721
       C 0.000000 -1.282877 -0.611721
       C 0.000000 0.000000 0.185210
       O 0.000000 0.000000 1.392088
       H 0.000000 2.133711 0.059851
       H -0.876575 1.319344 -1.256757
       H 0.876575 1.319344 -1.256757
       H 0.000000 -2.133711 0.059851
       H 0.876575 -1.319344 -1.256757
       H -0.876575 -1.319344 -1.256757
    &END COORD
    &TOPOLOGY
     &CENTER_COORDINATES T
     &END
    &END
    &KIND H
      BASIS_SET ORB DZVP-MOLOPT-PBE0-GTH-q1
      BASIS_SET AUX_FIT admm-dzp-q1
      POTENTIAL GTH-PBE0-q1
    &END KIND
    &KIND O
      BASIS_SET ORB DZVP-MOLOPT-PBE0-GTH-q6
      BASIS_SET AUX_FIT admm-dzp-q6
      POTENTIAL GTH-PBE0-q6
    &END KIND
    &KIND C
      BASIS_SET ORB DZVP-MOLOPT-PBE0-GTH-q4
      BASIS_SET AUX_FIT admm-dzp-q4
      POTENTIAL GTH-PBE0-q4
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

In the resulting output file, there is a TDDFPT section reporting the steps of the calculations.