Structure factor in Particle Mesh Ewald algorithm

1 September, 2021

This document aims to explain the ‘‘normalisation’’ of the structure factor $S_f$ in the Particle Mesh Ewald Algorithm. For background, one can consult this JCTC paper.

According to Coulomb’s law, the exact electrostatic energy of a charge distribution $\rho(\bf r)$, periodic under translations $\bf n$, is $$ E = \frac{1}{2} \iint \rho(\bf r) \sum_{\bf n} \frac{\rho ( \bf r’ + \bf n)}{|\bf r - \bf r^\prime|} \mathrm{d^3}{\bf r’} \mathrm{d^3}{\bf r} $$ Albeit its simple form, \ref{eq:Coulomb} is hardly practical for those who seek to capture long-ranged effects in organic macromolecules, where the integration is made unaffordable by infinite periodic images of a discontinuous distribution. Fortunately, long-ranged, periodic interaction calculations are bound in reciprocal space.

The charge distribution $\ket{\rho} $ is often better described by multipoles than by the real position basis ${\bra{\bf r}}$.

Spherical multipoles basis ${\bra{\delta_{l\mu}}}$ are a orthonormal basis in Hermitian space, as complete as any other, defined by $$ \begin{aligned}\delta_{\chi \mu } (\bf r - \bf R_a) &\equiv \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \delta (\bf r - \bf R_a)\\ \text{satisfying }\quad\bra{\delta_{l\mu}} \ket{\delta_{k\nu}} &= \delta_{lk} \delta_{\mu \nu}\end{aligned} $$ where $C_{l \mu}$ are the Spherical Tensor Gradient Operators (STGOs), and the ${}_a$ subscript denotes an atom. These basis provide a *spherical multipoles expansion* $$ \begin{aligned} q_{a, l\mu} &= \bra{\delta_{l\mu}} \ket{\rho_a} \\ \rho (\bf r) &= \sum_{a} \rho_a (\bf r - \bf R_a)\\ \rho (\bf r) &= \sum_{a} \bra{\bf r - \bf R_a} \sum_{l =0}^\infty \ket{\delta_{k\nu}} \bra{\delta_{l\mu}}\ket{\rho_a} \\ \rho (\bf r) &= \sum_{a} \sum_l^\infty q_{a, l \mu} \delta_{l \mu} ({\bf r - \bf R_a}) \end{aligned} $$ The higher multipole order, the more quickly Coulomb interaction decays away. Therefore, the charge distribution is often well represented by the lowest few orders of $l$. We extract a smooth part of the charge distribution to facilitate calculation in reciprocal space, by replacing the basis constructed from the Dirac delta with that from a normalised Gaussian $\chi$ of finite width $\frac{1}{\kappa}$: $$ \begin{aligned} \delta ( \bf r - \bf R_a) &\equiv \lim_{\kappa \to \infty} \chi (\bf r - \bf R_a)\\ \tilde{\rho} ( \bf r) &= \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \chi_{l \mu} ({\bf r - \bf R_a})\\ \rho (\bf r) &= \tilde{\rho} ( \bf r) + \Delta \rho(\bf r) \end{aligned} $$ where the correction part $\Delta \rho (\bf r)$ is left to be treated in real space. $$ \begin{gathered} \sum_{\bf n} \rho {\left( \bf r’ + \bf n \right)} = \underbrace{\sum_{\bf n} \tilde{\rho} ({ \bf r’ + \bf n})}_\text{smooth} + \underbrace{{\left[\sum_{\bf n} \rho {\left( \bf r’ + \bf n \right)} - \sum_{\bf n} \tilde{\rho} {\left( \bf r’ + \bf n \right)} \right]}}_ \text{discontinuous}\\ E = \frac{1}{2} \iint \rho(\bf r) \sum_{\bf n} \frac{\tilde{\rho} {\left( \bf r’ + \bf n \right)}}{{\left|\bf r - \bf r’ \right|}} \mathrm{d^3}{\bf r’} \mathrm{d^3}{\bf r} + \frac{1}{2} \iint \rho(\bf r) \sum_{\bf n} \frac{ \rho {\left( \bf r’ + \bf n \right)} - \tilde{\rho} {\left( \bf r’ + \bf n \right)}}{{\left|\bf r - \bf r’ \right|}} \mathrm{d^3}{\bf r’} \mathrm{d^3}{\bf r} \end{gathered} $$

\subsection{Baby steps toward reciprocal space energy} The Ewald summation is an reexpression of Coulomb’s law in the discrete wavevector space.

Periodicty is naturally incorporated into the smooth part charge density by spanning the gaussian with plane waves ${ \bf k}$, $$ \begin{gathered} \bf n = \sum_i^3 n_j \bf a_j\\ \bf A^*_i \cdot \bf a_j = 2 \pi \delta_{ij}\\ \bf k = \sum_i^3 m_i \bf A^*_i, m \in \mathbb{Z} \end{gathered} $$ where $\bf n$ is the aforementioned translational symmetry. $$ \begin{aligned} \tilde{\rho} ( \bf r) &= \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf r - \bf R_a} \ket{\chi } \\ \sum_{\bf n} \tilde{\rho} ( \bf r + \bf n) &= \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf r - \bf R_a} \frac{1}{V}\sum_{\bf k} \ket{\bf k} \bra{\bf k} \ket{\chi } \end{aligned} $$ After shifting some bra-kets around, the charge density can be reexpressed. $$ \begin{aligned} \sum_{\bf n} \tilde{\rho} ( \bf r + \bf n) &= \frac{1}{V}\sum_{\bf k} \bra{\bf k} \ket{\chi } \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \underbrace{\bra{\bf r - \bf R_a} \ket{\bf k}}_{\bra{\bf r} \ket{\bf k} \left(\bra{\bf R_a} \ket{\bf k} \right)^* } \\ \sum_{\bf n} \tilde{\rho} ( \bf r + \bf n) &= \frac{1}{V}\sum_{\bf k} \bra{\bf r} \ket{\bf k} \bra{\bf k} \ket{\chi } \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf k} \ket{\bf R_a} \\ \sum_{\bf n} \tilde{\rho} ( \bf r + \bf n) &= \frac{1}{V}\sum_{\bf k} \bra{\bf r} \ket{\bf k} \bra{\bf k} \ket{\chi } S_{\bf k} = \frac{1}{V}\sum_{\bf k} \bra{\bf k} \ket{\bf r}\bra{\chi }\ket{\bf k} S_{\bf k}^* \end{aligned} $$ where in the last line we used realness of $\rho$ and denoted $$ S_{\bf k} \equiv \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf k} \ket{\bf R_a} = \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (-i \bf k )}{(2l - 1)!!} \bra{\bf k} \ket{\bf R_a}. $$ as the structure factor. It is worth noticing $$ \sum_{\bf n} \rho( \bf r + \bf n) = \frac{1}{V}\sum_{\bf k} \bra{\bf k} \ket{\bf r} \left({ \lim_{\kappa \to \infty} \bra{\chi }\ket{\bf k}}\right) S_{\bf k}^* = \frac{1}{V}\sum_{\bf k} \bra{\bf k} \ket{\bf r} S_{\bf k}^* $$ because the Fourier transform of the Dirac delta is simply unity. $$ E_\text{recip} = \frac{1}{2} \iint \rho(\bf r) \sum_{\bf n} \frac{\tilde{\rho} \left({ \bf r’ + \bf n}\right)}{\left|{\bf r - \bf r’}\right|} \mathrm{d^3}{\bf r’} \mathrm{d^3}{\bf r} $$ Substituting \cref{eq:rhotil} and \cref{eq:rho} into reciprocal space part energy, we obtain $$ \begin{aligned} E_\text{recip}&= \frac{1}{2} \int \frac{1}{V}\sum_{\bf k} \bra{\bf k} \ket{\bf r} S_{\bf k}^* \frac{1}{V}\sum_{\bf k’} \bra{\bf k’} \ket{\chi } S_{\bf k’} \underbrace{\int \frac{\bra{\bf r’} \ket{\bf k’}}{{\bf r - \bf r’}} \mathrm{d^3}{\bf r’} }_{\bra{\bf r} \ket{\bf k’} \frac{4 \pi}{k^2}} \mathrm{d^3}{\bf r} \\ E_\text{recip}&= \sum_{\bf k, \bf k’} \frac{2\pi}{Vk'^2} \bra{\bf k} \underbrace{\frac{1}{V} \int \mathrm{d^3}{\bf r}\ket{\bf r} \bra{\bf r}}_1 \ket{\bf k’} S_{\bf k}^* \bra{\bf k’} \ket{\chi } S_{\bf k’}\\ E_\text{recip}&= \sum_{\bf k, \bf k’} \frac{2\pi}{Vk'^2} \delta_{\bf k’ \bf k} S_{\bf k}^* \bra{\bf k’} \ket{\chi } S_{\bf k’}\\ E_\text{recip}&= \boxed{\sum_{ \bf k} \frac{2\pi}{Vk^2} \bra{\bf k} \ket{\chi } { S_{\bf k}}^2} \end{aligned} $$ which is the Ewald expression of reciprocal space electrostatic energy. Given a desired relative error bound, only finite numbers of $\bf k$ need to be taken into account.

$$ \bf k = \sum_i^3 m_i \bf A^*_i,; m_i \in \left[\text{floor} - \frac{K_i - 1}{2}, \text{floor} \frac{K_i - 1}{2} \right] $$

The Ewald structure factor calculated in \cref{eq:esf} is computationally complex in that it scales as $O(N_a^2)$, where $N_a$ is the number of atoms in a unit cell of the system. The Particle Mesh Ewald method avoids this complexity by ‘‘meshing’’ the multipoles onto grid points. To make progress, we make all the following derivations, \emph{assuming $\bf m'$ is a set of meshgrid points with infinite resolution}, and then assess our conclusion in the sense of finite grid points $FFT$. $$ \begin{aligned} \sum_{\bf m’} \theta(\bf R_{\bf m’}) \xi \left({ \left({\bf r - \bf R_a}\right) - \bf R_{\bf m’}}\right) \approx \delta(\bf r -\bf R_a) \end{aligned} $$ Here, $\theta$ is Cardinal B-spline weighting function, normalised upon summation over grid points $m$ around the real atomic position $\bf R_a$, $\xi$ is the ``reciprocal’’ function of $\xi$ such that Dirac delta can be written as a convolution of $\xi$ and $\theta$. If the grid is resolved well enough, we can make $\bf R_{\bf m’} \to \bf R_{\bf m’} - \bf R_a$. Letting $\bf r = \bf R_{\bf m}$, $$ \begin{aligned} \sum_{\bf m’} \theta(\bf R_{\bf m’} - \bf R_a) \xi \left({\bf R_{\bf m} - \bf R_{\bf m’}}\right) = \delta(\bf R_{\bf m} - \bf R_a) \end{aligned} $$

Equipped with this expansion of Dirac delta, we can go back to \cref{eq:esf}. $$ \begin{aligned} S_{\bf k} &= \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf k} \sum_{\bf m} \ket{\bf R_{\bf m}} \bra{\bf R_{\bf m} - \bf R_a} \ket{ \delta}\\ S_{\bf k} &= \sum_{\bf m} \bra{\bf k} \ket{\bf R_{\bf m}} \sum_{\bf m’} \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \theta(\bf R_{\bf m’} - \bf R_a) \xi \left({\bf R_{\bf m} - \bf R_{\bf m’}}\right) \end{aligned} $$

Defining $$ Q(\bf R_{\bf m}) = \sum_{a} \sum_{l = 0}^{l_{\max}} q_{a, l \mu} \frac{C_{l \mu} (\nabla_a)}{(2l - 1)!!} \bra{\bf R_{\bf m} - \bf R_a} \ket{\theta} $$ We have $$ \begin{aligned} S_{\bf k} &= \sum_{\bf m} \bra{\bf k} \ket{\bf R_{\bf m}} \sum_{\bf m’} Q(\bf R_{\bf m’})\xi\left({\bf R_{\bf m} - \bf R_{\bf m’}}\right)\\ S_{\bf k} &= \sum_{\bf m} \bra{\bf k} \ket{\bf R_{\bf m}} {Q * \xi}\left({\bf R_{\bf m}}\right)\\ S_{\bf k} &= \mathcal{F} \left({{Q * \xi}\left({\bf R_{\bf m}}\right)}\right) =S^\text{PME}_{\bf k} \times \mathcal{F} \left({\xi\left({\bf R_{\bf m}}\right)}\right) \end{aligned} $$ where we have used convolution theorem and denoted $$ \begin{aligned} S^\text{PME}_{\bf k} &= \mathcal{F} \left({Q\left({\bf R_{\bf m}}\right)}\right). \end{aligned} $$ Noticing $$ \begin{gathered} \mathcal{F} \left({\xi\left({\bf R_{\bf m}}\right)}\right) \mathcal{F} \left({\theta\left({\bf R_{\bf m}}\right)}\right) = \mathcal{F} \left({\sum_{\bf m’} \theta(\bf R_{\bf m’}) \xi\left({\bf R_{\bf m} - \bf R_{\bf m’}}\right)}\right) = \mathcal{F} \left({\delta\left({\bf R_{\bf m}}\right)}\right) = 1\\ \mathcal{F} \left({\xi\left({\bf R_{\bf m}}\right)}\right) = \frac{1}{ \mathcal{F} \left({\theta\left({\bf R_{\bf m}}\right)}\right) } = \frac{1}{\theta_{\bf k}} \end{gathered} $$ We can conclude that the structure factor defined in \cref{eq:sfpme} differs from that in \cref{eq:esf} by a renormalisation factor $$ S_{\bf k} =\frac{S^\text{PME}_{\bf k} }{ \theta_{\bf k}} $$

$$ E^\text{PME}\text{recip}= \sum{ \bf k} \frac{2\pi}{Vk^2} \bra{\bf k} \ket{\chi } \left| \frac{S^\text{PME}{\bf k} }{ \theta{\bf k}}\right|^2 $$ \section{Formally deriving derivatives of the mesh function} The weight function $\theta$ is differentiated up to $l_\text{max}$ times. We need compact expressions for the evaluation of its derivatives.

$$ \theta = \prod_{d = 1}^3 M_n \left({ N_d ( \mathbf r - \mathbf R_a) \cdot \bf a^*_d + \frac{n}{2} }\right) $$ Denoting $u_d = N_d ( \mathbf r - \mathbf R_a) \cdot \bf a^*_d + \frac{n}{2} $ and $\partial_i \equiv \frac{\partial}{\partial u_i}$, we have for first derivative $$ \begin{aligned} \partial_i \theta &= \theta \partial_i \left({\sum_d \ln M_n(u_d) }\right)\\ &= \theta \sum_d \partial_i \ln M_n(u_d) \\ &= \theta \frac{M’_n(u_i)}{M_n(u_i)} \end{aligned} $$

where $M’_n, M’'_n$ are derivatives of $M_n$ derived in \texttt{derive_theta.py}. The above expression can be easily verified to be correct by writing it in the following form $$ \partial_i {(M_n(u_1) M_n(u_2) M_n(u_3))} = M_n(u_1) M_n(u_2) M_n(u_3) \frac{M’_n (u_i)}{M_n(u_i)}. $$ In numpy code, the fraction is a elementwise division.

Moving on to second derivative, $$ \begin{aligned} \partial_j \partial_i \theta &= \partial_j \left({\theta \partial_i \ln \theta}\right)\\ &= \frac{\left({\partial_i \theta}\right) \left({\partial_j \theta}\right)}{\theta} + \theta \partial_j \partial_i \ln M_n (u_i)\\ &= \frac{\left({\partial_i \theta}\right) \left({\partial_j \theta}\right)}{\theta} + \theta \delta_{ij} \left[{\frac{M’'_n (u_i)}{M_n(u_i)} - \left({\frac{M’_n(u_i)}{M_n(u_i)} }\right)^2}\right]\\ &= \theta \delta_{ij} \left[{\frac{M’'_n (u_i)}{M_n(u_i)} }\right] + \frac{\left({\partial_i \theta}\right) \left({\partial_j \theta}\right) - \delta_{ij} \left({\partial_i \theta}\right)^2 }{\theta} \end{aligned} $$

where the first term is diagonal and the second term is off-diagonal.

Derivations demonstrated above should allow us to express higher order derivatives of $\theta$ in terms of direct products of $M_d$ with more confidence and ease.

To get derivatives with respect to $\bf r$, chain rule can be employed $$ \frac{\partial }{\partial r_i} = \sum_j \frac{ \partial u_j}{\partial r_i} \partial_j = \sum_j - N_j A^*_{ji} \partial_j $$ where $A^*_{ji} \hat{\bf e}_i = \bf a^*_j$. Since $ - N_j A^*_{ji}$ are constant Jacobians, the transform is easy in arbitrary order of derivative.