Linear response with Casida and TDHF

Casida’s matrix formulation of linear response allows the calculation of the excitation energies of a finite system. For small molecules, this is normally the fastest way to obtain them. The very same electron-hole matrix problem describes both time-dependent density-functional theory (TDDFT) and time-dependent Hartree-Fock (TDHF); the latter is historically known as the random-phase approximation (RPA). In Octopus both are driven from CalculationMode=casida: whether the calculation is performed using TDDFT or TDHF is determined solely by the choice of functional or kernel (see theory level) of the underlying ground-state calculation.

Theory

The linear-response eigenvalue problem

Linear response recasts the neutral excitations of the interacting system as an eigenvalue problem in the space of single orbital transitions. We label each transition by a pair $ia$, meaning “occupied orbital $i$ $\to$ unoccupied (virtual) orbital $a$”, with Kohn-Sham (or Hartree-Fock) orbital energies $\varepsilon_i$, $\varepsilon_a$. The excitation energies $\omega_I$ and their eigenvectors $(\mathbf{X}_I,\mathbf{Y}_I)$ are the solutions of the non-Hermitian problem

$$ \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{A} \end{pmatrix} \begin{pmatrix} \mathbf{X}_I \\ \mathbf{Y}_I \end{pmatrix} = \omega_I \begin{pmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{pmatrix} \begin{pmatrix} \mathbf{X}_I \\ \mathbf{Y}_I \end{pmatrix}\,. $$

Here $\mathbf{X}$ collects the excitation ($i\to a$) and $\mathbf{Y}$ the de-excitation ($a\to i$) amplitudes. This is the standard polarization-propagator / RPA form; a self-contained pedagogic derivation (in which the method is called the RPA rather than TDHF) is given by Oddershede, Jørgensen and Yeager¹, while the density-functional version was derived to Casida²³. For the full-response and Tamm-Dancoff levels Octopus assumes real orbitals, so that $\mathbf{A}$ and $\mathbf{B}$ are real symmetric.

What distinguishes TDDFT from TDHF

The matrices $\mathbf{A}$ (the resonant block) and $\mathbf{B}$ (the coupling block) differ between the two theories only in the two-electron coupling. Writing the Coulomb integrals in Mulliken (charge-cloud) notation $(ia|jb)$, the TDHF (RPA) matrices are

$$ A_{ia,jb} = \delta_{ij}\delta_{ab}\,(\varepsilon_a - \varepsilon_i) + (ia|jb) - (ib|ja)\,, $$

$$ B_{ia,jb} = (ia|jb) - (ib|ja)\,, $$

where the negative terms are the non-local Fock exchange. The TDDFT (Casida) matrices instead replace that exchange by introducing the (semilocal) exchange-correlation kernel $f_{\rm xc}$ as,

$$ A_{ia,jb} = \delta_{ij}\delta_{ab}\,(\varepsilon_a - \varepsilon_i) + (ia|jb) + (ia|f_{\rm xc}|jb)\,, $$

$$ B_{ia,jb} = (ia|jb) + (ia|f_{\rm xc}|jb)\,. $$

This side-by-side comparison is exactly Eqns. 2 and 3 of Liang et al.⁴; hybrid functionals interpolate between the two by mixing a fraction of Fock exchange with $f_{\rm xc}$. In Octopus the TDHF matrices above are assembled for a ground state run at the Hartree-Fock theory level, and the TDDFT matrices for a DFT ground state.

The reduced Hermitian problem

For real orbitals and $(\mathbf{A}-\mathbf{B})$ positive definite, the $2N\times 2N$ non-Hermitian problem can be folded into an $N\times N$ Hermitian eigenvalue problem for $\omega_I^2$,

$$ (\mathbf{A}-\mathbf{B})^{1/2}(\mathbf{A}+\mathbf{B})(\mathbf{A}-\mathbf{B})^{1/2}\,\mathbf{F}_I = \omega_I^2\,\mathbf{F}_I\,, \qquad \mathbf{F}_I \propto (\mathbf{A}-\mathbf{B})^{-1/2}(\mathbf{X}_I+\mathbf{Y}_I)\,, $$

which is the form derived by Stratmann, Scuseria and Frisch⁵ and the one solved by Octopus for the full TDHF (RPA) problem. A loss of positive-definiteness of $(\mathbf{A}-\mathbf{B})$ signals a (triplet) instability of the reference state, in which case only the Tamm-Dancoff level is available, cf. the next subsection.

In pure (semilocal) TDDFT the exchange-correlation contributions to $\mathbf{A}$ and $\mathbf{B}$ are identical, so $(\mathbf{A}-\mathbf{B})$ reduces to the diagonal matrix of orbital-energy differences and the reduced problem takes Casida’s celebrated form²

$$ \Omega_{ia,jb} = \delta_{ij}\delta_{ab}\,(\varepsilon_a-\varepsilon_i)^2 + 2\sqrt{\varepsilon_a-\varepsilon_i}\;K_{ia,jb}\;\sqrt{\varepsilon_b-\varepsilon_j}\,, \qquad \mathbf{\Omega}\,\mathbf{F}_I = \omega_I^2\,\mathbf{F}_I\,, $$

with the coupling kernel $K_{ia,jb} = (ia|jb) + (ia|f_{\rm xc}|jb)$.

Approximations and oscillator strengths

Neglecting the coupling block ($\mathbf{B}=\mathbf{0}$) yields the Tamm-Dancoff approximation (TDA), $\mathbf{A},\mathbf{X}_I = \omega_I,\mathbf{X}_I$ ⁶; applied on top of a Hartree-Fock reference this is the configuration-interaction-singles (CIS) method. Keeping only the diagonal (or degenerate blocks) of $\mathbf{A}$ gives the single-pole Petersilka approximation⁷, and the bare orbital-energy differences $\omega_{ia}=\varepsilon_a-\varepsilon_i$ correspond to the independent-particle level.

For each excitation the dipole oscillator strength, from which the absorption spectrum is built, is

$$ f_I = \frac{2}{3}\,\omega_I \sum_{n \in \{x,y,z\}} \left| \langle 0 | \hat{r}_n | I \rangle \right|^2\,, $$

where $|0\rangle$ is the ground state, $\hat{r}$ is the dipole transition operator, and the excited state $|I\ket$ is reconstructed from the eigenvectors $\mathbf{X}_I$ and $\mathbf{Y}_I$.

Theory levels in Octopus

The CasidaTheoryLevel variable selects which of these levels to compute (several may be requested at once, since they share most of the matrix construction):

For a Hartree-Fock ground state only the tamm_dancoff and full lrtddft_casida (i.e. TDHF/RPA) levels are available; the petersilka and variational levels are skipped, and the calculation is currently restricted to the non-distributed matrix layout, spin-unpolarized/collinear spin, and no photon modes. Hybrid functionals are not yet supported.

Running a calculation

To perform a Casida (or TDHF) calculation you first need a ground-state calculation, and a calculation of unoccupied states; then you run the code with CalculationMode=casida.

A tutorial is available here: tutorial on Casida calculations.

Examples of Casida with photons, or for magnetic clusters can be found at share/testsuite/linear_response/07-casida-photons.* and share/testsuite/linear_response/09-casida-gga.* . A worked TDHF (RPA) example is provided at share/testsuite/linear_response/12-casida_hf.* .

References