# Tutorial on the smooth particle-mesh Ewald algorithm

#### Abstract.

These notes explain the smooth particle-mesh Ewald algorithm (SPME) in detail. The accompanying library libpme6 is a fully-featured implementation of SPME that can be used for experimentation and reference. This content can also be found in PDF format here. Suggestions and feedback will be appreciated.

## 1. Introduction

### 1.1. Motivation and overview of the SPME algorithm

Consider two atoms, labeled 1 and 2, contained in a parallelepiped with sides determined by the (positively oriented) vectors $\mathbf{\ell}_{1},\mathbf{\ell}_{2},\mathbf{\ell}_{3}\in\mathbb{R}^{3}$ (see Figure 1).

Then the matrix $\mathtt{L}=[\mathbf{\ell}_{1},\mathbf{\ell}_{2},\mathbf{\ell}_{3}]\in\mathbb{R% }^{3\times 3}$ is such that $V=\det\mathtt{L}>0$. Now suppose that atom $i$ is at a position $\mathbf{r}_{i}=\mathtt{L}\,\mathbf{x}_{i}$ where $\mathbf{x}_{i}\in[0,1)^{3}$ for $i=1,2$ and that $r_{ij}=\|\mathbf{r}_{i}-\mathbf{r}_{j}\|$ (see Figure 2).

Between these two atoms there is a force that depends on their distance and whose corresponding potential energy, $U_{ij}=U_{ij}(r_{ij})$, is shown in Figure 3.

If instead of just two atoms we had a total of $N>2$ atoms (see Figure 4), then the total potential energy would be equal to

 $H=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}U_{ij}(r_{ij}).$ (1.1)

By solving Newton’s second law11Force equals mass times acceleration. we would be able to simulate the molecular system of $N$ atoms that interact according to the (pairwise) forces given by $F=-\nabla U$. In practice, whatever value of $N$ that can be tractable with current computers is still far away from the more realistic figure of $N_{A}=6.022\times 10^{23}$ (Avogadro’s number) which would capture the real physics of the system.

One trick to somehow get closer to $N_{A}$ is to use a tractable value of $N$ and impose periodic boundary conditions (PBCs). This means that we now have an infinite lattice of simulation boxes with identical copies of the particles in each box. By computing the interactions between the atoms in the original box and their copies (only the copies in the neighboring simulation boxes are typically taken into account —this is called the minimum image convention—), we can make our simulations more realistic.

The total potential energy (1.1) then becomes

 $H=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}U_{ij}(\|\mathbf{% r}_{i}-\mathbf{r}_{j}\|)+\frac{1}{2}\sum_{\stackrel{\mathbf{n}\in\mathbb{Z}^{3% }}{\mathbf{n}\neq{0}}}\sum_{i=1}^{N}\sum_{j=1}^{N}U_{ij}(\|\mathbf{r}_{i}-% \mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|).$

The naive approach to PBCs requires $\mathcal{O}(N^{2})$ calculations and is unfeasible, hence the interest of finding efficient ways of carrying out this task. Ewald summation [EWA21] takes advantage of ideas from Harmonic analysis to reduce the complexity to $\mathcal{O}(N^{3/2})$ operations, which is still impractical if we want to solve the equations of motion of a large system. The Particle-Mesh Ewald (PME) method [DYP93] uses the Fast Fourier Transform to bring the complexity down to $\mathcal{O}(N\log N)$, and this algorithm has been a staple of Molecular Simulation since its appearance. The Smooth Particle-Mesh Ewald (SPME) improves on PME by giving a sufficiently smooth energy function whose derivative can be obtained analytically, which results in more realistic simulations due to improved energy conservation.

Some of the people involved in the development of the SPME algorithm appear in Figure 6.

## 2. Smooth Particle-Mesh Ewald

We derive the formulas for Ewald summation following [WIL06] and then we introduce the smooth Particle-Mesh Ewald method from [EPBDLP95].

### 2.1. Decomposition of the potential energy function

Ewald summation can be applied whenever the pairwise potential energy $U=U(r)$ is proportional to $r^{-n}$ for some $n\in\mathbb{N}$ (a concrete instance of $U$ appears in Figure 7). In that case, we can write $U$ as a sum of the form

 $U(r)=U_{0}(r)+U_{1}(r),$ (2.1)

where the short-range term, $U_{0}$, decays fast as $r$ goes to infinity (so that it can be truncated) and the long-range term, $U_{1}$, is very spread out and smooth, so that its Fourier transform is concentrated (and can also be truncated).

There are many ways of splitting the potential energy but [ND57] establishes the optimality for a class of functions involving the (incomplete) Gamma function. If we define

 $\varphi_{n}(r)=1-\frac{\Gamma(n/2,\alpha^{2}\,\pi\,r^{2})}{\Gamma(n/2)},$

where $\alpha^{2}$ is a parameter to be discussed later, then we clearly have

 $\frac{1}{r^{n}}=\frac{1}{r^{n}}\,\left(\varphi_{n}(r)+(1-\varphi_{n}(r))\right).$

We are going to concern ourselves with the case $n=6$ (see Figure 7), and to define $\psi(r)=r^{-6}\varphi_{6}(r)$.

The free parameter $\alpha$ is determined by numerically solving the equation

 $r^{-6}\left(1-\varphi_{6}(r_{\text{cutoff}};\alpha)\right)=\varepsilon$

with prescribed values of $r_{\text{cutoff}}$ and $\varepsilon$ (the function $1-\varphi_{n}$ must be decreasing for $r>r_{\text{cutoff}})$. This makes the short-range part of the computation involving effectively vanish for $r>r_{\text{cutoff}}$.

Recall that in the case of $N>2$ particles, we have

 $H=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}U_{ij}(\|\mathbf{% r}_{i}-\mathbf{r}_{j}\|)+\frac{1}{2}\sum_{\stackrel{\mathbf{n}\in\mathbb{Z}^{3% }}{\mathbf{n}\neq{0}}}\sum_{i=1}^{N}\sum_{j=1}^{N}U_{ij}(\|\mathbf{r}_{i}-% \mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|).$

where $r_{ij}$ is the distance between atoms $i$ and $j$. We can split the pairwise potential $U_{ij}$ above as in (2.1).

### 2.2. Ewald summation

The summation of the long-range terms (with $b_{i}>0$),

 $U_{1}=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}b_{i}b_{j}\,% \,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|)+\frac{1}{2}\sum_{\stackrel{\mathbf{n% }\in\mathbb{Z}^{3}}{\mathbf{n}\neq{0}}}\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}% \,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|)$ (2.2)

can be rearranged as

 $U_{1}=\frac{1}{2}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\sum_{i=1}^{N}\sum_{j=1}^{N% }b_{i}b_{j}\,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|)-% \frac{1}{2}\sum_{i=1}^{N}b_{i}^{2}\,\psi(0),$ (2.3)

where $\psi(0)=\lim_{r\to 0^{+}}\psi(r)$.

#### Lemma 2.1 (Poisson’s summation formula).

If $D(\mathbf{r})=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\delta(\mathbf{r}-\mathtt{L}% \mathbf{n})$, then $\widehat{D}(\mathbf{h})=\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3}}\delta(% \mathbf{h}-\mathtt{L}^{-\mathtt{T}}\mathbf{d})$.

#### Lemma 2.2.

Consider the distributions $\rho(\mathbf{r})=\sum_{i=1}^{N}b_{i}\,\delta(\mathbf{r}-\mathbf{r}_{i})$, $\eta(\mathbf{r})=\rho(-\mathbf{r})$, and $D(\mathbf{r})=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\delta(\mathbf{r}-\mathtt{L}% \mathbf{n})$. The following identity holds:

 $\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\,\psi% (\|\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|)=\int\psi(r)\,\left(% \rho\star\eta\star D\right)(\mathbf{r})\,\mathrm{d}\mathbf{r},$

where $r=\|\mathbf{r}\|$.

#### Proof.

First,

 $\displaystyle\rho\star\eta\,(\mathbf{r})$ $\displaystyle=\int\rho(\mathbf{s})\,\eta(\mathbf{r}-\mathbf{s})\,\mathrm{d}% \mathbf{s}$ $\displaystyle=\int\rho(\mathbf{s})\,\rho(\mathbf{s}-\mathbf{r})\,\mathrm{d}% \mathbf{s}$ $\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\int\delta(\mathbf{s}-% \mathbf{r}_{i})\,\delta(\mathbf{r}-\mathbf{s}-\mathbf{r}_{j})\,\mathrm{d}% \mathbf{s}$ $\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\,\delta(\mathbf{r}-(% \mathbf{r}_{i}-\mathbf{r}_{j})).$

By the associative property of the convolution, we have

 $\rho\star\eta\star D=\left(\rho\star\eta\right)\star D,$

so

 $\displaystyle\rho\star\eta\star D\,(\mathbf{r})$ $\displaystyle=\int\rho\star\eta\,(\mathbf{s})\,D(\mathbf{s}-\mathbf{r})\,% \mathrm{d}\mathbf{s}$ $\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\int\delta(\mathbf{s}-(% \mathbf{r}_{i}-\mathbf{r}_{j}))\,D(\mathbf{s}-\mathbf{r})\,\mathrm{d}\mathbf{s}$ $\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\,D(\mathbf{s}-(\mathbf{r}% _{i}-\mathbf{r}_{j}))$ $\displaystyle=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\sum_{i=1}^{N}\sum_{j=1}^{N}b_% {i}b_{j}\,\delta(\mathbf{s}-(\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n% })).$

Finally,

 $\displaystyle\int\psi(r)\,\left(\rho\star\eta\star D\right)(\mathbf{r})\,% \mathrm{d}\mathbf{r}$ $\displaystyle=\int\psi(r)\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\sum_{i=1}^{N}\sum_% {j=1}^{N}b_{i}b_{j}\,\delta(\mathbf{s}-(\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{% L}\mathbf{n}))\,\mathrm{d}\mathbf{r}$ $\displaystyle=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\sum_{i=1}^{N}\sum_{j=1}^{N}b_% {i}b_{j}\,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|).$

At this point we can go back to (2.3) and rewrite it as

 $U_{1}=\frac{1}{2}\int\psi(r)\,\left(\rho\star\eta\star D\right)(\mathbf{r})\,% \mathrm{d}\mathbf{r}-\frac{1}{2}\sum_{i=1}^{N}b_{i}^{2}\,\psi(0),$ (2.4)

Furthermore, using Plancherel’s Theorem [WOL03] and the fact that $\rho\star\eta\star D$ is real-valued, we can write

 $\int\psi(r)\,\left(\rho\star\eta\star D\right)(\mathbf{r})\,\mathrm{d}\mathbf{% r}=\int\widehat{\psi}(h)\,\widehat{\left(\rho\star\eta\star D\right)}(\mathbf{% h})\,\mathrm{d}\mathbf{h},$ (2.5)

where $h=\|\mathbf{h}\|$. Consequently, equation (2.4) becomes

 $U_{1}=\frac{1}{2}\int\widehat{\psi}(h)\,\widehat{\left(\rho\star\eta\star D% \right)}(\mathbf{h})\,\mathrm{d}\mathbf{h}-\frac{1}{2}\sum_{i=1}^{N}b_{i}^{2}% \,\psi(0),$ (2.6)

By the properties of the Fourier transform of a convolution [WOL03] and Poisson’s summation formula 2.1,

 $\widehat{\left(\rho\star\eta\star D\right)}(\mathbf{h})=\widehat{\rho}(\mathbf% {h})\,\widehat{\rho}(-\mathbf{h})\,\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3% }}\delta(\mathbf{h}-\mathtt{L}^{-\mathtt{T}}\mathbf{d})=\big|\widehat{\rho}(% \mathbf{h})\big|^{2}\,\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3}}\delta(% \mathbf{h}-\mathtt{L}^{-\mathtt{T}}\mathbf{d}).$

From the above and (2.5), we have

 $\displaystyle\int\widehat{\psi}\,\widehat{\left(\rho\star\eta\star D\right)}$ $\displaystyle=\int\widehat{\psi}(h)\,\left(\big|\widehat{\rho}(\mathbf{h})\big% |^{2}\,\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3}}\delta(\mathbf{h}-\mathtt{% L}^{-\mathtt{T}}\mathbf{d})\right)\,\mathrm{d}\mathbf{h}$ $\displaystyle=\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3}}\widehat{\psi}(\|% \mathtt{L}^{-\mathtt{T}}\mathbf{d}\|)\,\big|\widehat{\rho}(\mathtt{L}^{-% \mathtt{T}}\mathbf{d})\big|^{2}$ $\displaystyle=\frac{1}{V}\sum_{\mathbf{m}}\widehat{\psi}(m)\,\big|\widehat{% \rho}(\mathbf{m})\big|^{2}$ $\displaystyle=\frac{1}{V}\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\big|% \widehat{\rho}(\mathbf{m})\big|^{2}+\frac{1}{V}\lim_{m\to 0}\widehat{\psi}(m)% \,\big|\widehat{\rho}({0})\big|^{2}$ $\displaystyle=\frac{1}{V}\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\big|% \widehat{\rho}(\mathbf{m})\big|^{2}+\frac{1}{V}\lim_{m\to 0}\widehat{\psi}(m)% \,\left(\sum_{i=1}^{N}b_{i}\right)^{2},$

where $\mathbf{m}=\mathtt{L}^{-\mathtt{T}}\mathbf{d}$ and $m=\|\mathbf{m}\|$. So, to sum up, we have arrived at the identity:

 $\displaystyle U_{1}$ $\displaystyle=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}b_{i}% b_{j}\,\,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|)+\frac{1}{2}\sum_{\stackrel{% \mathbf{n}\in\mathbb{Z}^{3}}{\mathbf{n}\neq{0}}}\sum_{i=1}^{N}\sum_{j=1}^{N}b_% {i}b_{j}\,\psi(\|\mathbf{r}_{i}-\mathbf{r}_{j}+\mathtt{L}\mathbf{n}\|)$ $\displaystyle=\frac{1}{2V}\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\big|% \widehat{\rho}(\mathbf{m})\big|^{2}+\frac{1}{2V}\widehat{\psi}(0)\,\left(\sum_% {i=1}^{N}b_{i}\right)^{2}-\frac{1}{2}\psi(0)\,\sum_{i=1}^{N}b_{i}^{2},$

### 2.3. Fast Fourier Transform and Particle-Mesh method

Ewald summation using optimal parameters requires $\mathcal{O}(N^{3/2})$ operations [dPS80] but it can be modified so that it involves only about $\mathcal{O}(N\log N)$ operations by using the Fast Fourier Transform [EPBDLP95, HE88].

We now concern ourselves with summing (and later taking derivatives of) the expression

 $\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\big|\widehat{\rho}(\mathbf{m})\big% |^{2}.$ (2.7)

By the definition of the (inverse) Fourier Transform, we see that

 $\widehat{\rho}(\mathbf{m})=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,% \mathbf{m}^{\mathtt{T}}\mathbf{r}_{j}},$

where $\mathbf{r}_{j}=\mathtt{L}\,\mathbf{x}_{j}$ for some $\mathbf{x}_{j}\in[0,1)^{3}$. Moreover,

 $\widehat{\rho}(\mathbf{m})=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,% \mathbf{m}^{\mathtt{T}}\mathbf{r}_{j}}=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi% \mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathtt{L}^{-1}\,\mathtt{L}\,\mathbf{x}_{j}% }=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}% \mathbf{x}_{j}}.$

We are going to approximate the complex exponential terms above by means of Cardinal B-Splines [SCH73]. An order $p$ B-spline function is the $p$-fold convolution of the uniform density in the unit interval (see Figure 9).

By construction, each order $p$ B-spline is of class $C^{p-1}$, their integral equals 1, and they have compact support.

#### Remark.

Observe that these splines play the role of approximations to the identity when used in a meshing of the unit cell. It is in this sense that they play a role in approximating the Dirac delta functions in $\rho$.

Now, to discretize the simulation box we consider $K_{1},K_{2},K_{3}\in\mathbb{N}$ and let

 $\mathtt{K}=\left[\begin{smallmatrix}K_{1}&&\\ &K_{2}&\\ &&K_{3}\end{smallmatrix}\right].$

For each $\mathbf{r}=\mathtt{L}\mathbf{x}$, we set $\mathbf{u}=\mathtt{K}\mathbf{x}$ and we have

 $\displaystyle\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{m}^{\mathtt{T}}\mathbf{r}}$ $\displaystyle=\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathbf{x}}=% \mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathtt{K}^{-1}\mathbf{u}}% \approx\sum_{\boldsymbol{\mathbf{\ell}}\in\mathbb{Z}^{3}}W(\mathbf{u}-% \boldsymbol{\mathbf{\ell}})\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T% }}\mathtt{K}^{-1}\boldsymbol{\mathbf{\ell}}}$ $\displaystyle=\sum_{\mathbf{k}}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}W(\mathbf{u}-% (\mathbf{k}+\mathtt{K}\mathbf{n}))\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{% \mathtt{T}}\mathtt{K}^{-1}(\mathbf{k}+\mathtt{K}\mathbf{n})}$ $\displaystyle=\sum_{\mathbf{k}}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}W(\mathbf{u}-% (\mathbf{k}+\mathtt{K}\mathbf{n}))\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{% \mathtt{T}}\mathtt{K}^{-1}\mathbf{k}}\underbrace{\mathrm{e}^{2\pi\mathrm{i}\,% \mathbf{d}^{\mathtt{T}}\mathbf{n}}}_{=1}$ $\displaystyle=\sum_{\mathbf{k}}\sum_{\mathbf{n}\in\mathbb{Z}^{3}}W(\mathbf{u}-% (\mathbf{k}+\mathtt{K}\mathbf{n}))\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{% \mathtt{T}}\mathtt{K}^{-1}\mathbf{k}},$

where we have used the Euclidean algorithm to write

 $\boldsymbol{\mathbf{\ell}}=\mathbf{k}+\mathtt{K}\mathbf{n}\quad\text{with}% \quad k_{i}\in\left\{0,1,\ldots,K_{i}-1\right\}\quad\text{and}\quad\mathbf{n}% \in\mathbb{Z}^{3}.$

We also define the forward (unnormalized) DFT of the 3D array $Q_{\mathbf{k}}$ as

 $\widehat{Q}_{\mathbf{d}}=\text{DFT}[Q_{\mathbf{k}}]_{\mathbf{d}}=\sum_{k_{1}=0% }^{K_{1}-1}\sum_{k_{2}=0}^{K_{2}-1}\sum_{k_{3}=0}^{K_{3}-1}Q_{\mathbf{k}}\,% \mathrm{e}^{-2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathtt{K}^{-1}\mathbf{k}}.$

And, likewise, the backward DFT is

 $\text{IDFT}[\widehat{Q}_{\mathbf{d}}]_{\mathbf{k}}=\sum_{d_{1}=0}^{K_{1}-1}% \sum_{d_{2}=0}^{K_{2}-1}\sum_{d_{3}=0}^{K_{3}-1}\widehat{Q}_{\mathbf{d}}\,% \mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathtt{K}^{-1}\mathbf{k}}.$

Using the previous results, we write

 $\displaystyle\widehat{\rho}(\mathbf{m})$ $\displaystyle=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{% \mathtt{T}}\mathbf{x}_{j}}$ $\displaystyle\approx\sum_{\mathbf{k}}\underbrace{\sum_{j=1}^{N}b_{j}\sum_{% \mathbf{n}\in\mathbb{Z}^{3}}W(\mathbf{u}-(\mathbf{k}+\mathtt{K}\mathbf{n}))}_{% =Q_{\mathbf{k}}}\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{\mathtt{T}}\mathtt{K}^% {-1}\mathbf{k}}$ $\displaystyle=\sum_{\mathbf{k}}Q_{\mathbf{k}}\,\mathrm{e}^{2\pi\mathrm{i}\,% \mathbf{d}^{\mathtt{T}}\mathtt{K}^{-1}\mathbf{k}}=\widehat{Q}_{\mathbf{d}}.$

#### Remark.

It is important to notice that since the splines have compact support, the infinite series above is actually a finite sum.

#### Remark.

The following diagram depicts the relationships between the different objects involved in the SPME method:

 $\begin{CD}\rho @>{\approx}>{}>Q\\ @V{\text{IFT}}V{}V@V{}V{\text{IDFT}}V\\ \widehat{\rho}@>{\approx}>{}>\widehat{Q}\end{CD}$

Going back to (2.7), we carry out the following approximation:

 $\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\big|\widehat{\rho}(\mathbf{m})\big% |^{2}=\sum_{\mathbf{m}\neq{0}}\widehat{\psi}(m)\,\widehat{\rho}(-\mathbf{m})\,% \widehat{\rho}(\mathbf{m})\approx\sum_{\mathbf{d}\neq{0}}\widehat{\psi}(\|% \mathtt{L}^{-\mathtt{T}}\,\mathbf{d}\|)\,\widehat{Q}_{-\mathbf{d}}\,\widehat{Q% }_{\mathbf{d}}$

Now define

 $\mathbf{h}\neq 0$

with

 $0\leq d_{i}<\frac{1}{2}K_{i}$

(note that $-\frac{K_{i}}{2}\leq h_{i}<\frac{K_{i}}{2}-1$). Then,

 $\sum_{d_{1}=-\frac{K_{1}}{2}}^{\frac{K_{1}}{2}-1}\sum_{d_{2}=-\frac{K_{2}}{2}}% ^{\frac{K_{2}}{2}-1}\sum_{d_{3}=-\frac{K_{3}}{2}}^{\frac{K_{3}}{2}-1}\widehat{% \psi}(\|\mathtt{L}^{-\mathtt{T}}\,\mathbf{d}\|)\,\widehat{Q}_{-\mathbf{d}}\,% \widehat{Q}_{\mathbf{d}}=\sum_{d_{1}=0}^{K_{1}-1}\sum_{d_{2}=0}^{K_{2}-1}\sum_% {d_{3}=0}^{K_{3}-1}A_{\mathbf{d}}\,\widehat{Q}_{-\mathbf{d}}\,\widehat{Q}_{% \mathbf{d}}.$

Applying the following property of the Fourier transform (see [WOL03]),

 $\int\widehat{f}\,\mathrm{d}g=\int\widehat{g}\,\mathrm{d}f,$

we write (using the notation of the previous section and keeping in mind that the result works in the distributional sense too)

 $\int{\widehat{\psi}\,\widehat{\eta}}\,\widehat{\rho}=\int\widehat{\widehat{% \psi}\,\widehat{\eta}}\,\rho=\int\left(\widehat{C}\star\rho\right)\,\rho$

where we have used that $\widehat{fg}=\widehat{f}\star\widehat{g}$ and that, for $f=f(x)$, $\widehat{\widehat{f}}(x)=f(-x)$ (again, see [WOL03]).

Therefore, the algorithm boils down to:

 $\displaystyle Q_{\mathbf{k}}$ $\displaystyle\rightarrow\widehat{Q}_{-\mathbf{d}}=\text{DFT}(Q)_{d}$ $\displaystyle\rightarrow A_{\mathbf{d}}\cdot\text{DFT}(Q)_{\mathbf{d}}$ $\displaystyle\rightarrow\text{IDFT}(A\cdot\text{DFT}(Q))_{\mathbf{k}}=(\text{% IDFT}(A)\star Q)_{\mathbf{k}}$ $\displaystyle\rightarrow(\text{IDFT}(A)\star Q)_{\mathbf{k}}\cdot Q_{\mathbf{k% }},$

## Appendix A Assorted results in Harmonic analysis

#### Theorem A.3 ([WOL03]).

For absolutely continuous measures $\mu$, $\nu$ in $\mathbb{R}^{n}$. We have

 $\int\hat{\mu}\,\mathrm{d}\nu=\int\hat{\nu}\,\mathrm{d}\mu$

## References

• [DYP93] T. Darden, D. York and L. Pedersen (1993) Particle mesh Ewald: An $N\log(N)$ method for Ewald sums in large systems. The Journal of Chemical Physics 98 (12), pp. 10089. External Links: Document, ISSN 00219606, Link Cited by: 1.1.
• [dPS80] S. W. de Leeuw, J. W. Perram and E. R. Smith (1980-10) Simulation of Electrostatic Systems in Periodic Boundary Conditions. I. Lattice Sums and Dielectric Constants. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 373 (1752), pp. 27–56. External Links: Document, ISSN 1364-5021, Link Cited by: 2.3.
• [EPBDLP95] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen (1995) A smooth particle mesh Ewald method. The Journal of Chemical Physics 103 (19), pp. 8577. External Links: Document, ISSN 00219606, Link Cited by: 2.3, 2.
• [EWA21] P. P. Ewald (1921) Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik 369 (3), pp. 253–287. External Links: Document, ISSN 1521-3889, Link Cited by: 1.1.
• [HE88] R. W. Hockney and J. J. W. Eastwood (1988) Computer Simulation Using Particles. Adam Hilger. External Links: ISBN 9780852743928, Link Cited by: 2.3.
• [ND57] B. R. A. Nijboer and F. W. De Wette (1957) On the calculation of lattice sums. Physica 23 (1-5), pp. 309–321. External Links: Document, ISSN 00318914, Link Cited by: 2.1.
• [SCH73] I. J. Schoenberg (1973) Cardinal Spline Interpolation. Society for Industrial and Applied Mathematics. External Links: Document, ISBN 978-0-89871-009-0, Link Cited by: 2.3.
• [WIL06] D. E. Williams (2006) Accelerated convergence treatment of ${R}^{-n}$ lattice sums. in U. Shmueli (Ed.), International Tables for Crystallography, Vol. B, pp. 385–397. Cited by: 2.
• [WOL03] T. H. Wolff (2003) Lectures on harmonic analysis. University Lecture Series, Vol. 29, American Mathematical Society, Providence, RI. Note: With a foreword by Charles Fefferman and preface by Izabella Laba, Edited by Izabella Laba and Carol Shubin External Links: ISBN 0-8218-3449-5, Link Cited by: A.3, 2.2, 2.3.