deletions | additions
diff --git a/Casida, Tamm-Dancoff, and excited-state forces.tex b/Casida, Tamm-Dancoff, and excited-state forces.tex
index c9ae23e..8480c5e 100644
--- a/Casida, Tamm-Dancoff, and excited-state forces.tex
+++ b/Casida, Tamm-Dancoff, and excited-state forces.tex
...
\end{align}
and the absorption cross-section is
\begin{align}
\sigma_{ij} \left( \omega \right) = \frac{4 \pi \omega}{c} \tilde{\alpha} \ \mathcal{I} \alpha_{ij} \left( \omega
\right) \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 \sigma} 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 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} \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} 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| v + f_{\rm xc} \left| \varphi_c \right> \left| \varphi_v
\right> \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| v + f_{\rm xc} \left| \varphi_c \right> \left| \varphi_v
\right> \right>\ .
\end{align} (Octopus currently only supports $f_{\rm xc}$ for LDA-type functionals.)
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}.
...
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} \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
...
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} 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} x_{cv}\,.
\end{align}
The dipole matrix elements are
\begin{align}
\vec{d}_k = \sum_{cv} x_{cv}\, \vec{d}_{cv} / \sqrt{\epsilon_c -
\epsilon_c} \epsilon_c}\ .
\end{align}
An alternate approach for finding excitation energies is to look for many-body eigenstates
...
\left[ \begin{array}{c|c}
A & B \\ \hline
B^{*} & -A
\end{array} \right] \vec{x} = \omega
\vec{x} \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} \vec{x}\ .
\end{align}
The Tamm-Dancoff approximation to these equations is identical to the ordinary TDDFT Tamm-Dancoff approximation.
Note that all
of the levels of theory we have discussed use the same Coulomb and $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
...
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) \right)\ , \\
\varphi ^{\rm T} = \frac{1}{\sqrt{2}} \left( \varphi_{c \uparrow} \varphi_{v \uparrow}
- \varphi_{c \downarrow} \varphi_{v \downarrow}
\right) \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}. excitations~\cite{Onida2002}
\begin{align}
\left< \varphi^{\rm S} \right|
v \hat{v} +
f_{\rm \hat{f}_{\rm xc} \left| \varphi^{\rm S} \right>
= \left< \varphi \right|
v \hat{v} +
f_{\rm \hat{f}_{\rm xc}^{\uparrow \uparrow} +
f_{\rm \hat{f}_{\rm xc}^{\uparrow \downarrow} \left| \varphi \right> = \left< \varphi \right|
v \hat{v} + 2 f_{\rm xc} \left| \varphi \right> \\
\left< \varphi^{\rm T} \right|
v \hat{v} +
f_{\rm \hat{f}_{\rm xc} \left| \varphi^{\rm T} \right>
= \left< \varphi \right|
f_{\rm \hat{f}_{\rm xc}^{\uparrow \uparrow} -
f_{\rm \hat{f}_{\rm xc}^{\uparrow \downarrow} \left| \varphi
\right> \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.
...
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(
r \vec{r} \right) \varphi_{v'} \left(
r \vec{r} \right) P \left[ \varphi_c \varphi_v
\right] dr \right]\ .
\end{align}
Our basic parallelization strategy for computation of the matrix elements is by domains, as discussed in section~\ref{sec:parallelization},
...
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}. 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 again is 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.
...
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 as in previous 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. energy
\begin{align}
F^{k}_{\partial R_{i \alpha}} = -\frac{\partial E_k}{\partial R_{i \alpha}} = F_{\partial R_{i \alpha}} - \frac{\partial \omega_k}{\partial R_{i
\alpha}} \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{align}
\frac{\partial \omega_k}{\partial R_{i \alpha}}
= &= \left< x_k \left| \frac{\partial A}{\partial R_{i \alpha}} \right| x_k \right> \\
\left<\varphi_c \varphi_v \left| \frac{\partial A}{\partial R_{i \alpha}} \right| \varphi_{c'} \varphi_{v'} \right>
= &=
\left<\varphi_c \left| \frac{\partial H}{\partial R_{i \alpha}} \right| \varphi_{c'} \right> \delta_{vv'} \\ \nonumber
- \left<\varphi_v \left| \frac{\partial H}{\partial R_{i \alpha}} \right| \varphi_{v'} \right> \delta_{cc'}
+ \left<\varphi_c \varphi_v \left| K_{\rm xc} \frac{\partial \rho}{\partial R_{i \alpha}} \right| \varphi_{c'} \varphi_{v'} \right>
\end{align}
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,
for leading to a significant savings in computational time
compared to the calculation of the excited states.
