\section{The Sternheimer formulation of linear-response}
\label{sec:sternheimer}
%David, Xavier

In textbooks, perturbation theory is formulated in terms of sums over
states and response functions. These are useful theoretical
constructions that permit a good description and understanding of the
underlying physics. However, this is not always a good starting point for
numerical applications, since it involves the calculation of a large
number of eigenvectors, infinite sums over these eigenvectors, and functions that depend on two or more spatial variables.

An interesting approach that avoids the problems
mentioned above is the formulation of perturbation theory
in terms of differential equations for the variation of the
wave-functions. In the literature, this is usually called
the Sternheimer equation~\cite{Sternheimer_1951} or density functional perturbation theory (DFPT)~\cite{Baroni_2001}. Although a
perturbative technique, it avoids the use of empty states, and has a
favorable scaling with the number of atoms.

Octopus implements a generalized version of the Sternheimer equation that is
able to cope with both static and dynamic response in and out of
resonance~\cite{Andrade_2007}. The method is suited for linear and
non-linear response; higher-order Sternheimer equations can be
obtained for higher-order variations. For second-order response,
however, we apply the \(2n\,+\,1\) theorem (also known as Wigner's  \(2n\,+\,1\) rule)~\cite{Gonze_1989,Corso_1996} to get the
coefficients directly from first-order response variations.

In the Sternheimer formalism, we consider the response to 
a monochromatic perturbative field \(\lambda
\delta{\hat v}(\vec{r})\cos\left(\omega{t}\right)\). This perturbation induces a variation in the time-dependent KS orbitals, which we denote $\delta\varphi_{n}(\vec r, \omega)$. These variations allow us to calculate response observables, for example, the frequency-dependent polarizability. 

In order to calculate the variations of the orbitals we need to solve a linear equation that only depends on the occupied orbitals (atomic units are used throughout)
%\begin{multline}
\begin{equation}
  \label{eq:sternheimer}
  \left\{\hat H  - \epsilon_n\pm\omega +
    \mathrm{i}\eta\right\}\delta\varphi_{n}(\vec r, \pm\omega) = -\mathrm{\hat P}_c\,\delta{\hat H}(\pm\omega) \varphi_n(\vec r)
  \,,
\end{equation}
%
where the variation of the time-dependent density, given by
%
%\begin{multline}
\begin{equation}
  \label{eq:varrho}
  \delta{n}(\vec r, \omega) = \sum_k f_k \Big\{
    \left[\varphi_n(\vec{r})\right]^*\delta\varphi_{n}(\vec r, \omega)\\
     + \left[\delta\varphi_{n}(\vec r, -\omega)\right]^*\varphi_n(\vec{r})
  \Big\}\ ,
\end{equation}
%\end{multline}
needs to be calculated self-consistently. The first-order variation of the KS Hamiltonian is
\begin{multline}
  \label{eq:h1}
  \delta{\hat H}(\omega)=
  \delta{\hat v}(\vec{r})
  +\int \mathrm{d}\vec{r}' \frac{\delta{n}(\vec{r}',\omega)}{|\vec{r}-\vec{r}'|}\\
  +\int \mathrm{d}\vec{r}' f_{\rm xc}(\vec r, \vec r', \omega)\,\delta{n}(\vec{r}', \omega)
  \ .
\end{multline}
%
\(\mathrm{\hat P}_c\) is a projection operator, and
\(\eta\) a positive in\-fi\-ni\-te\-si\-mal, essential to obtain the
correct position of the poles of the causal response function, and
consequently to obtain the imaginary part of the
po\-la\-ri\-za\-bi\-li\-ty and remove the divergences of
  the equation for resonant frequencies. In the usual implementation of DFPT, $\hat P_c = 1 - \sum^{\rm occ}_n \left| \varphi_n \right> \left< \varphi_n \right|$ which effectively removes the components of \(\delta\varphi_{n}(\vec r, \pm\omega)\) in the
subspace of the occupied ground-state wave-functions. In linear
response, these components do not contribute to the variation of the
density.

We have found that it is not strictly necessary to project out the occupied subspace, 
the crucial part is simply to remove the
projection of $\delta \varphi_n$ on $\varphi_n$ (and any other states degenerate with it),
which is not physically meaningful and arises from a phase convention. To fix the phase, it is
sufficient to apply a minimal projector
$\hat P_n = 1 - \sum_m^{\epsilon_m = \epsilon_n} \left| \varphi_m \right> \left< \varphi_m \right|$.
We optionally use this approach to obtain the entire response wavefunction, not just the projection in the
unoccupied subspace, which is needed for obtaining effective masses in $\vec{k} \cdot \vec{p}$ theory. While the full projection can become time-consuming
for large systems, it saves time overall since it increases the condition number of the matrix for the linear solver,
and thus reduces the number of solver iterations required to attain a given precision.

We also have implemented the Sternheimer formalism when non-integer occupations are used, as appropriate for metallic systems. In this case weighted projectors are added to both sides of eq.~(\ref{eq:sternheimer})~\cite{DeGironcoli}. We have generalized the equations to the dynamic case \cite{Strubbethesis}. The modified Sternheimer equation is
\begin{multline}
  \left\{\hat H  - \epsilon_n\pm\omega +
    \mathrm{i}\eta + \sum_m \alpha_m \left| \varphi_m \right> \left< \varphi_m \right| \right\}\delta\varphi_{n}(\vec r, \pm\omega) = \\ 
    -\left[ \tilde{\theta}_{{\rm F},n} - \sum_m \beta_{n,m} \left| \varphi_m \right> \left< \varphi_m \right| \right] \delta{\hat H}(\pm\omega) \varphi_n(\vec r)
  \,,
\end{multline}
where
%
\begin{align}
\alpha_n &= {\rm max} \left( \epsilon_{\rm F} - 3 \sigma - \epsilon_n, 0 \right)\ , \\
\beta_{n,m} &= \tilde{\theta}_{{\rm F},n} \tilde{\theta}_{n,m} + \tilde{\theta}_{{\rm F},m} \tilde{\theta}_{m,n} + \alpha_m \frac{\tilde{\theta}_{{\rm F},n} - \tilde{\theta}_{n,m}}{\epsilon_n - \epsilon_m \mp \omega} \tilde{\theta}_{m,n}\ ,
\end{align}
%
$\sigma$ is the broadening energy, and $\tilde{\theta}_{ij}$ is the smearing scheme's approximation to the Heaviside function $\theta \left( \left( \epsilon_i - \epsilon_j \right) / \sigma \right)$.
Apart from semiconducting smearing (\textit{i.e.} the original equation above, which corresponds to the zero temperature limit),
the code offers the standard Fermi-Dirac~\cite{Mermin_1965}, Methfessel-Paxton~\cite{Methfessel_1989}, spline~\cite{Holender_1995}, and cold~\cite{Marzari_1999} smearing schemes.
Additionally, we have developed a scheme for handling arbitrary fractional occupations, which do not have to be defined by a function of the energy eigenvalues \cite{Strubbethesis}.

In order to solve eq.~(\ref{eq:sternheimer}) we use a self-consistent iteration scheme similar to the one used for ground-state DFT. In each iteration we need to solve a sparse linear problem where the operator to invert is the shifted KS Hamiltonian. For real wavefunctions and a real shift (as for the static case), we can use conjugate gradients. When the shift is complex, a non-Hermitian iterative solver is required. We have found that a robust and efficient solver is the quasi-minimal residual (QMR) method~\cite{Freund_1991}.
