This alternative formulation of the forces can be extended to obtain the second-order derivatives of the energy with respect to the atomic displacements~\cite{Andrade2010thesis}, which are required to calculate vibrational properties as discussed in section~\ref{sec:sternheimer}. In general, the perturbation operator associated with an ionic displacement can be written as

\begin{equation}
  \label{eq:ionicpertmod}
  \frac{\partial {v}_\alpha(\vec{r}-
    \vec{R}_\alpha)}{\partial R_{i\alpha}} = 
  -{v}_{\alpha}(\vec{r}-
    \vec{R}_\alpha)\frac{\partial}{\partial r_i}
  -\frac{\partial}{\partial
    r_i}{v}_\alpha(\vec{r}-
    \vec{R}_\alpha)\ .
\end{equation}
%
Using this expression, the terms of the dynamical matrix, eq.~(\ref{eq:dynmatrix}), are evaluated as
%
\begin{multline}
  \left< \varphi_n \left| \frac{\partial \hat{v}_{\alpha}}{\partial R_{i\alpha}} \right| \frac{\partial \varphi_n}{\partial R_{j \beta}} \right> = -\left[
 \left< \varphi_n \left| \hat{v}_{\alpha} \right| \frac{\partial^2 \varphi_n}{\partial R_{j \beta} \partial r_i} \right> \right.\\ + \left. \left< \frac{\partial \varphi_n}{\partial r_i}\left| \hat{v}_{\alpha} \right| \frac{\partial \varphi_n}{\partial R_{j \beta}} \right>\right] + {\rm c.c.}\ ,
\end{multline}
and
\begin{multline} 
  \left< \varphi_n \left| \frac{\partial^2 \hat{v}_{\alpha}}{\partial R_{i\alpha} \partial R_{j\alpha}} \right| \varphi_n \right> =
\left[  \left< \frac{\partial^2 \varphi_n}{\partial r_i \partial r_j} \left| \hat{v}_{\alpha} \right| \varphi_n \right>\right.\\ + 
 \left. \left< \frac{\partial \varphi_n}{\partial r_i} \left| \hat{v}_{\alpha}\right| \frac{\partial \varphi_n}{\partial r_j} \right>\right] + {\rm c.c.}\ .
\end{multline}
