\section{Linear response in the electron-hole basis}
\label{sec:casida}

An alternate approach to linear response is not to solve for the response function but rather for its poles (the excitation energies $\omega_k$) and residues (\textit{e.g.} electric dipole matrix elements $\vec{d}_k$) \cite{Strubbe2012}. The polarizability is given by
\begin{align}
  \alpha_{ij} \left( \omega \right) = \sum_k \left[ \frac{\left( \hat{i} \cdot \vec{d}_k \right)^{*} \left( \hat{j} \cdot \vec{d}_k \right)}{\omega_k - \omega - i \delta}
    + \frac{\left( \hat{i} \cdot \vec{d}_k \right)^{*} \left( \hat{j} \cdot \vec{d}_k \right)}{\omega_k + \omega + i \delta} \right]
\end{align}
and the absorption cross-section is
\begin{align}
\sigma_{ij} \left( \omega \right) = \frac{4 \pi \omega}{c} \tilde{\alpha} \  \mathrm{Im}\ \alpha_{ij} \left( \omega \right)\ ,
\end{align}
where $\tilde{\alpha}$ is the fine-structure constant.
%and the oscillator strengths are
%\begin{align}
%f_k = \frac{2 m \omega_k}{\hbar^2} \left| d_k \right|^2
%\end{align}
The simplest approximation to use is the random-phase approximation (RPA), in which the excitation energies are given by the differences of unoccupied and occupied KS eigenvalues, $\omega_{c v} = \epsilon_c - \epsilon_v$. The corresponding dipole matrix elements are $\vec{d}_{cv} = \left< \varphi_c \left| \vec{r} \right| \varphi_v \right>$ \cite{Onida2002}. (As implemented in the code, this section will refer only to the case of a system without partially occupied levels.)

The RPA is not a very satisfactory approximation, however. The full solution within TDDFT is given by a non-Hermitian matrix eigenvalue equation, with a basis consisting of both occupied-unoccupied ($v \rightarrow c$) and unoccupied-occupied ($c \rightarrow v$) KS transitions. The equation reads as
\begin{align}
  \left[ \begin{array}{c|c}
               A      & B \\ \hline
               -B^{*} & -A
    \end{array} \right] \vec{x} = \omega \vec{x}\ ,
\end{align}
where the $A$ matrices couple $v \rightarrow c$ transitions among themselves and $c \rightarrow v$ among themselves, while the $B$ matrices couple the two types of transitions. They have the form~\cite{Onida2002}
\begin{align}
  \left< \varphi_{c'} \right| \left< \varphi_{v'} \right| A \left| \varphi_c \right> \left| \varphi_v \right> &=
  \left( \epsilon_c - \epsilon_v \right) \delta_{cc'} \delta_{vv'} \\ \nonumber
  &+ \left< \varphi_{c'} \right| \left< \varphi_{v'} \right| \hat{v}_{\rm c} + \hat{f}_{\rm xc} \left| \varphi_c \right> \left| \varphi_v \right>\ , \\
  \left< \varphi_{c'} \right| \left< \varphi_{v'} \right| B \left| \varphi_c \right> \left| \varphi_v \right> &=
  \left< \varphi_{c'} \right| \left< \varphi_{v'} \right| \hat{v}_{\rm c} + \hat{f}_{\rm xc} \left| \varphi_c \right> \left| \varphi_v \right>\ .
\end{align}
where $\hat{v}_{\rm c}$ is the Coulomb kernel, and $\hat{f}_{\rm xc}$ is the exchange-correlation kernel. Currently, only the kernel for LDA-type functionals is supported in Octopus. The implementation for more advanced functionals is planned for future releases.

We do not solve the full equation in Octopus, but provide a hierarchy of approximations. An example calculation for the N$_2$ molecule with each theory level is shown in Table \ref{tab:nitrogen_casida}.
The lowest approximation we use is RPA.  The next is the single-pole approximation of Petersilka \textit{et al.} \cite{Petersilka1996},
in which only the diagonal elements of the matrix are considered.
Like in the RPA case, the eigenvectors and dipole matrix elements are simply the KS transitions. The positive eigenvalues are $\omega_{cv} = \epsilon_c - \epsilon_v + A_{cvcv}$.
This can be a reasonable approximation when there is little mixing between KS transitions,
but generally fails when there are degenerate or nearly degenerate transitions.

A next level of approximation is the Tamm-Dancoff approximation to TDDFT \cite{Hirata1999}
in which the B blocks are neglected and thus we need only consider the
occupied-unoccupied transitions. The matrix equation is reduced to a Hermitian problem
of half the size of the full problem:
\begin{align}
A \vec{x} = \omega \vec{x}\ .
\end{align}
Interestingly, the Tamm-Dancoff approximation is often found to give superior results to
the full solution, for example for molecular potential-energy surfaces
%cite one of the forces papers? e.g. Tsukagoshi, Sitt, Hutter etc.
or when hybrid functionals are used, which can suffer from a ``triplet instability''
in which the lowest triplet state is lower in energy than the ground state \cite{Casida2009}.
The dipole matrix elements are now a superposition of the KS ones:
\begin{align}
\vec{d}_k = \sum_{cv} \vec{d}_{cv}\, x_{cv}\ .
\end{align}

When the wavefunctions are real, the full problem can be collapsed into a Hermitian one of the same
size as the Tamm-Dancoff matrix, known as Casida's equation \cite{Casida1995,Jamorski1996}.
\begin{align}
  \left( \epsilon_c - \epsilon_v \right)^2 \delta_{cc'} \delta_{vv'} + 2 \sqrt{\epsilon_{c'} - \epsilon_{v'}}
  B_{cvc'v'} \sqrt{\epsilon_c - \epsilon_v} = \omega^2 x_{cv}\,.
\end{align}
The dipole matrix elements are
\begin{align}
\vec{d}_k = \sum_{cv} \vec{d}_{cv}\, x_{cv}\, \sqrt{\frac{\epsilon_c - \epsilon_v}{\omega_k}}\ .
\end{align}

An alternate approach for finding excitation energies is to look for many-body eigenstates
of the DFT Hamiltonian which are orthogonal to the ground state. In the ``second-order constrained
variational'' or CV(2) theory \cite{Ziegler2009},
second-order perturbation theory from the ground-state density
yields equations quite similar to the linear-response approach, despite their different origin:
\begin{align}
  \left[ \begin{array}{c|c}
               A      & B \\ \hline
               B^{*} & -A
    \end{array} \right] \vec{x} = \omega \vec{x}\ .
\end{align}
We implement the case of real wavefunctions and eigenvectors, in which case (as for Casida's equation) a Hermitian matrix equation
for only the occupied-unoccupied transitions can be written:
\begin{align}
\left( A + B \right) \vec{x} = \omega \vec{x}\ .
\end{align}
The Tamm-Dancoff approximation to these equations is identical to the ordinary TDDFT Tamm-Dancoff approximation.

Note that all the levels of theory we have discussed use the same Coulomb and $\hat{f}_{\rm xc}$ matrix elements, so the code can calculate the results for multiple levels of theory with a small extra effort. We can also consider alternative perturbations in this framework beyond the dipole approximation for properties such as inelastic X-ray scattering \cite{Sakko2010}.

\begin{table}
\centering
\begin{tabular}{ccccc|c}
RPA & Petersilka & TDA & Casida & CV(2) & Exp't \\ \hline
8.234    & 9.421      &  9.343       & 9.254  & 9.671       & 9.309 \\
8.234    & 9.421      &  9.343       & 9.254  & 10.279      & 9.309 \\
9.671    & 9.671      &  9.671       & 9.671  & 10.279      & 9.921 \\
9.671    & 10.241     &  10.237      & 10.221 & 10.792      & 10.270 \\
9.671    & 10.245     &  10.241      & 10.224 & 10.801      & 10.270 \\
9.671    & 11.028     &  10.931      & 10.921 & 11.077      & 12.199 \\
\end{tabular}
\caption{
\label{tab:nitrogen_casida}
The first 6 excitation energies (in eV) for the N$_2$ molecule with different approximations to TDDFT in the electron-hole basis: the random phase approximation (RPA), Petersilka, Tamm-Dancoff approximation (TDA), Casida and CV(2). The VWN LDA parametrization~\cite{Vosko_1980} was used for the exchange-correlation functional, the bond length is 1.098 \AA, the real-space grid was a sphere of radius 7.4 \AA~with spacing 0.16 \AA, and 16 unoccupied states were used. The experimental data is from Ref.~\cite{Ben-Shlomo_1990}.}
\end{table}

For a non-spin-polarized system, the excitations separate into a singlet and a triplet subspace,
which are superpositions of singlet and triplet KS transitions:
\begin{align}
  \varphi ^{\rm S} = \frac{1}{\sqrt{2}} \left( \varphi_{c \uparrow} \varphi_{v \uparrow}
  + \varphi_{c \downarrow} \varphi_{v \downarrow} \right)\ , \\
  \varphi ^{\rm T} = \frac{1}{\sqrt{2}} \left( \varphi_{c \uparrow} \varphi_{v \uparrow}
  - \varphi_{c \downarrow} \varphi_{v \downarrow} \right)\ .
\end{align}
The signs are reversed from the situation for a simple pair of electrons, since we are instead dealing
with an electron and a hole. There are of course two other triplet excitations ($m = \pm 1$) which are degenerate with the $m = 0$ one above.
Rather than performing spin-polarized ground-state and linear-response calculations, we can use
the symmetry between the spins in a non-spin-polarized system to derive a form of the kernel to use
in obtaining singlet and triplet excitations~\cite{Onida2002}
\begin{align}
  \left< \varphi^{\rm S} \right| \hat{v}_{\rm c} + \hat{f}_{\rm xc} \left| \varphi^{\rm S} \right>
  = \left< \varphi \right| \hat{v}_{\rm c} + \hat{f}_{\rm xc}^{\uparrow \uparrow} + \hat{f}_{\rm xc}^{\uparrow \downarrow} \left| \varphi \right> = \left< \varphi \right| \hat{v}_{\rm c} + 2 \hat{f}_{\rm xc} \left| \varphi \right> \\
  \left< \varphi^{\rm T} \right| \hat{v}_{\rm c} + \hat{f}_{\rm xc} \left| \varphi^{\rm T} \right>
  = \left< \varphi \right| \hat{f}_{\rm xc}^{\uparrow \uparrow} - \hat{f}_{\rm xc}^{\uparrow \downarrow} \left| \varphi \right>\ .
\end{align}
These kernels can be used in any of the levels of theory above: RPA, Petersilka, Tamm-Dancoff, Casida, and CV(2).
The corresponding electric dipole matrix elements are as in the spin-polarized case for singlet excitations.
For triplet excitations, they are identically zero, and only higher-order electromagnetic processes can excite them.

There are three main steps in the calculation: calculation of the matrix, diagonalization of the matrix, and calculation
of the dipole matrix elements. The first step generally takes almost all the computation time, and is the most important
to optimize. Within that step, the Coulomb part (since it is non-local) is much more time-consuming than the $\hat{f}_{\rm xc}$
part. We calculate it by solving the Poisson equation (as for the Hartree potential) for each column of the matrix,
to obtain a potential P for the density $\varphi_c \left( r \right)^{*} \varphi_v \left( r \right)$, and then for each row
computing the matrix element as
\begin{align}
  \left< \varphi_{c'} \varphi_{v'} \right| v \left| \varphi_c \varphi_v \right>
  = \int \mathrm{d}\vec{r} \varphi_{c'} \left( \vec{r} \right) \varphi_{v'} \left( \vec{r} \right) P \left[ \varphi_c \varphi_v \right]\ .
\end{align}

Our basic parallelization strategy for computation of the matrix elements is by domains, as discussed in section~\ref{sec:parallelization},
but we add an additional level of
parallelization here over occupied-unoccupied pairs. We distribute the columns of the matrix, and do not distribute
the rows, to avoid duplication of Poisson solves. We can reduce the number of matrix elements to be computed by
almost half using the Hermitian nature of the matrix, \textit{i.e.} $M_{cv,c'v'} = M_{c'v',cv}^{*}$. If there are $N$
occupied-unoccupied pairs, there are $N$ diagonal matrix elements, and the $N \left( N - 1 \right) / 2$ remaining
off-diagonal matrix elements are distributed as evenly as possible among the columns. If $N-1$ is even,
there are $\left( N - 1 \right) / 2$ for each column; if $N-1$ is odd, half of the columns have $N/2 - 1$ and
half have $N/2$. See Fig. \ref{fig:casida_parallelization} for examples of the distribution. The columns then are assigned to the available
processors in a round-robin fashion. The diagonalization step is performed by direct diagonalization with LAPACK~\cite{LAPACK}
in serial; since it generally accounts for only a small part of the computation time, parallelization of this step is not very important. The final step is calculation of
the dipole matrix elements, which amounts to only a small part of the computation time, and uses only domain parallelization.
Note that the triplet kernel lacks the Coulomb term, and so is considerably faster to compute.

Using the result of a calculation of excited states by one of these methods, and a previous calculation
of vibrational modes with the Sternheimer equation, we can compute forces in each excited state, which can be used 
for excited-state structural relaxation or molecular dynamics \cite{Strubbe_forces}. Our formulation allows us to do this without introducing any extra summations over empty states, unlike previous force implementations \cite{Sitt2007,Tsukagoshi2012,Hutter2003}.
The energy of a given excited state $k$ is a sum of the ground-state energy and the excitation energy: $E_k = E_0 + \omega_k$.
The force is then given by the ground-state force, minus the derivative of the excitation energy:
\begin{align}
  F^{k}_{i \alpha} = -\frac{\partial E_k}{\partial R_{i \alpha}} = F_{i \alpha} - \frac{\partial \omega_k}{\partial R_{i \alpha}}\ .
\end{align}
Using the Hellman-Feynman Theorem we find the last term without introducing any additional sums over unoccupied states.
In the particular case of the Tamm-Dancoff approximation we have
\begin{equation}
  \frac{\partial \omega_k}{\partial R_{i \alpha}} = \left< x_k \left| \frac{\partial \hat{A}}{\partial R_{i \alpha}} \right| x_k \right>\ ,
  \end{equation}
  and
  \begin{multline}
  \left<\varphi_c \varphi_v \left| \frac{\partial \hat{A}}{\partial R_{i \alpha}} \right| \varphi_{c'} \varphi_{v'} \right> =
  \left<\varphi_c \left| \frac{\partial \hat{H}}{\partial R_{i \alpha}} \right| \varphi_{c'} \right> \delta_{vv'}  \\
  - \left<\varphi_v \left| \frac{\partial \hat{H}}{\partial R_{i \alpha}} \right| \varphi_{v'} \right> \delta_{cc'}
  + \left<\varphi_c \varphi_v \left| \hat{K}_{\rm xc} \frac{\partial \rho}{\partial R_{i \alpha}} \right| \varphi_{c'} \varphi_{v'} \right>\,.
\end{multline}
Analogous equations apply for the difference of eigenvalues, Petersilka, and CV(2) theory levels. (The slightly more complicated Casida case has not yet been implemented.)
The Coulomb term, with no explicit dependence on the atomic positions, does not appear, leading to a significant savings in computational time
compared to the calculation of the excited states.