deletions | additions       diff --git a/Sternheimer.tex b/Sternheimer.tex  index 5c8b7b8..983f047 100644  --- a/Sternheimer.tex  +++ b/Sternheimer.tex 
...
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 can apply the \(2n\,+\,1\) theorem (also known as Wigner's \(2n\,+\,1\) rule)~\cite{Gonze_1989,Baroni_2001,Corso_1996} 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   monochromatic perturbative field \(\lambda  \delta{V}(\vec{r})\cos\left(\omega{t}\right)\). This perturbation induces a variation in the time-dependent Kohn-Sham orbitals, that we denote $\delta\varphi_{k}(\vec r, \omega)$. These variations allows to calculate response quantities like the frequency-dependent polarization. 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} 
...
\Big\}\ ,  \end{equation}  %\end{multline}  need needs  to be calculated self-consistently.%\end{multline}  The first order variation of the Kohn-Sham Hamiltonian is %\begin{multline}  \begin{equation}  \label{eq:h1} 
...
\delta{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{equation}  %  Where \(\mathrm{P}_c\) is the a  projector onto the unoccupied subspace 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 obtain the imaginary part of the  po\-la\-ri\-za\-bi\-li\-ty and remove the divergences of  this the  equation for resonant frequencies. The projector \(\mathrm{P}_c\) effectively removes In  the usual implementation of DFTP, $P_c = 1 - \sum^{\rm occ}_k \left| \varphi_k \right> \left< \varphi_k \right|$ which effectively removesthe  components of \(\delta\varphi_{k}(\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, and therefore we can safely ignore the  projector for first-order response calculation. This is important for  large systems as the cost of the calculation of the projections scales  quadratically with the number of orbitals. density.  In fact We have found that  it is not strictly necessary to project out the occupied subspace with  $P_c = 1 - \sum^{\rm occ}_k \left| \varphi_k \right> \left< \varphi_k \right|$. The subspace,   the  crucial part is simply to remove the projection of $\delta \varphi_k$ on $\varphi_k$ (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  $P_k = 1 - \sum_l^{\epsilon_l = \epsilon_k} \left| \varphi_l \right> \left< \varphi_l \right|$.  not valid for smearing.  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. 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, 
...
We also have implemented the Sternheimer equation with smearing, as appropriate for metallic systems,  in which weighted projectors are added to both sides of the equation. %cite De Geroncoli.  % FIXME eqn here  Apart from semiconducting (\textit{i.e.} the original equation above, which is the limit T $\rightarrow$ 0), zero temperature limit),  the code offers the standard Fermi-Dirac, Methfessel-Paxton, Fermi-Dirac~\cite{Mermin_1965}, Methfessel-Paxton~\cite{Methfessel_1989},  spline, and cold 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 my thesis, and D. A. Strubbe and S. G. Louie, in prep. 
...
%M Methfessel and AT Paxton, Phys. Rev. B 40, 3616 (1989).  %JM Holender, MJ Gillan, MC Payne, and AD Simpson, Phys. Rev. B 52, 967 (1995).  In order to solve equation~\ref{eq:sternheimer} we use a self-consistent iteration scheme similar to the ones used for ground-state DFT. Each iteration we need to solve a linear problem where the operator to invert is the shifted Kohn-Sham Hamiltonian, as the shift is complex, a non-Hermitian solver is required. This operator is sparse, so an iterative solver is required. We have found that a robust solver is the quasi-minimal residual (QMR) method~\cite{Freund_1991}.  We can solve for linear response to various different perturbations. To compute vibrational frequencies,  we calculate response to the ionic displacement perturbation  $\partial H/\partial R_{i \alpha} = \partial V_{\alpha}/\partial x_i$, for each direction $i$ and atom $\alpha$.