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 1,2,33subscriptnormal-ℓ1subscriptnormal-ℓ2subscriptnormal-ℓ3superscript3\mathbf{\ell}_{1},\mathbf{\ell}_{2},\mathbf{\ell}_{3}\in\mathbb{R}^{3} (see Figure 1).

Figure 1. Simulation box spanned by three vectors.

Then the matrix 𝙻=[1,2,3]3×3𝙻subscriptnormal-ℓ1subscriptnormal-ℓ2subscriptnormal-ℓ3superscript33\mathtt{L}=[\mathbf{\ell}_{1},\mathbf{\ell}_{2},\mathbf{\ell}_{3}]\in\mathbb{R% }^{3\times 3} is such that V=det𝙻>0V𝙻0V=\det\mathtt{L}>0. Now suppose that atom iii is at a position 𝐫i=𝙻𝐱isubscript𝐫i𝙻subscript𝐱i\mathbf{r}_{i}=\mathtt{L}\,\mathbf{x}_{i} where 𝐱i[0,1)3subscript𝐱isuperscript013\mathbf{x}_{i}\in[0,1)^{3} for i=1,2i12i=1,2 and that rij=𝐫i-𝐫jsubscriptrijnormsubscript𝐫isubscript𝐫jr_{ij}=\|\mathbf{r}_{i}-\mathbf{r}_{j}\| (see Figure 2).

Figure 2. Two atoms at a distance rijsubscriptrijr_{ij}.

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

Figure 3. Decomposition of the Lennard Jones 12-6 potential into repulsive (powers of -1212-12) and dispersive (powers of -66-6) contributions.

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

H=12i=1Nijj=1NUij(rij).H12superscriptsubscripti1Nsuperscriptsubscriptsuperscriptijj1NsubscriptUijsubscriptrijH=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N}U_{ij}(r_{ij}). (1.1)
Figure 4. Argon gas (N=103Nsuperscript103N=10^{3} atoms).

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

Figure 5. The protein DHFR solvated in water (N=23536N23536N=23536 atoms, so N2108superscriptN2superscript108N^{2}\approx 10^{8} and Nlog2N105Nsubscript2Nsuperscript105N\log_{2}N\approx 10^{5}).

One trick to somehow get closer to NAsubscriptNAN_{A} is to use a tractable value of NNN 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=12i=1Nijj=1NUij(𝐫i-𝐫j)+12𝐧0𝐧3i=1Nj=1NUij(𝐫i-𝐫j+𝙻𝐧).H12superscriptsubscripti1Nsuperscriptsubscriptsuperscriptijj1NsubscriptUijnormsubscript𝐫isubscript𝐫j12subscriptsuperscript𝐧0𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1NsubscriptUijnormsubscript𝐫isubscript𝐫j𝙻𝐧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 𝒪(N2)𝒪superscriptN2\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 𝒪(N3/2)𝒪superscriptN32\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 𝒪(NlogN)𝒪NN\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.

(a) Paul Peter Ewald (1888-1985).
(b) Tom Darden and Lee Pedersen.
(c) Mike Crowley.
Figure 6. Some of the people behind SPME. Pictures reproduced from http://journals.iucr.org/a/issues/2012/01/00/wx0006/ and http://www.psc.edu/science/Darden/.

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)UUrU=U(r) is proportional to r-nsuperscriptrnr^{-n} for some nnn\in\mathbb{N} (a concrete instance of UUU appears in Figure 7). In that case, we can write UUU as a sum of the form

U(r)=U0(r)+U1(r),UrsubscriptU0rsubscriptU1rU(r)=U_{0}(r)+U_{1}(r), (2.1)

where the short-range term, U0subscriptU0U_{0}, decays fast as rrr goes to infinity (so that it can be truncated) and the long-range term, U1subscriptU1U_{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

φn(r)=1-Γ(n/2,α2πr2)Γ(n/2),subscriptφnr1normal-Γn2superscriptα2πsuperscriptr2normal-Γn2\varphi_{n}(r)=1-\frac{\Gamma(n/2,\alpha^{2}\,\pi\,r^{2})}{\Gamma(n/2)},

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

1rn=1rn(φn(r)+(1-φn(r))).1superscriptrn1superscriptrnsubscriptφnr1subscriptφnr\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=6n6n=6 (see Figure 7), and to define ψ(r)=r-6φ6(r)ψrsuperscriptr6subscriptφ6r\psi(r)=r^{-6}\varphi_{6}(r).

Figure 7. Decomposition of the function r-6superscriptr6r^{-6} into short-range and long-range terms.

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

r-6(1-φ6(rcutoff;α))=εsuperscriptr61subscriptφ6subscriptrcutoffαεr^{-6}\left(1-\varphi_{6}(r_{\text{cutoff}};\alpha)\right)=\varepsilon

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

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

H=12i=1Nijj=1NUij(𝐫i-𝐫j)+12𝐧0𝐧3i=1Nj=1NUij(𝐫i-𝐫j+𝙻𝐧).H12superscriptsubscripti1Nsuperscriptsubscriptsuperscriptijj1NsubscriptUijnormsubscript𝐫isubscript𝐫j12subscriptsuperscript𝐧0𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1NsubscriptUijnormsubscript𝐫isubscript𝐫j𝙻𝐧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 rijsubscriptrijr_{ij} is the distance between atoms iii and jjj. We can split the pairwise potential UijsubscriptUijU_{ij} above as in (2.1).

2.2. Ewald summation

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

U1=12i=1Nijj=1Nbibjψ(𝐫i-𝐫j)+12𝐧0𝐧3i=1Nj=1Nbibjψ(𝐫i-𝐫j+𝙻𝐧)subscriptU112superscriptsubscripti1Nsuperscriptsubscriptsuperscriptijj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j12subscriptsuperscript𝐧0𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j𝙻𝐧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

U1=12𝐧3i=1Nj=1Nbibjψ(𝐫i-𝐫j+𝙻𝐧)-12i=1Nbi2ψ(0),subscriptU112subscript𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j𝙻𝐧12superscriptsubscripti1Nsuperscriptsubscriptbi2ψ0U_{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 ψ(0)=limr0+ψ(r)ψ0subscriptnormal-→rsuperscript0ψr\psi(0)=\lim_{r\to 0^{+}}\psi(r).

Lemma 2.1 (Poisson’s summation formula).

If D(𝐫)=𝐧3δ(𝐫-𝙻𝐧)D𝐫subscript𝐧superscript3δ𝐫𝙻𝐧D(\mathbf{r})=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\delta(\mathbf{r}-\mathtt{L}% \mathbf{n}), then D^(𝐡)=1V𝐝3δ(𝐡-𝙻-𝚃𝐝)normal-^D𝐡1Vsubscript𝐝superscript3δ𝐡superscript𝙻𝚃𝐝\widehat{D}(\mathbf{h})=\frac{1}{V}\sum_{\mathbf{d}\in\mathbb{Z}^{3}}\delta(% \mathbf{h}-\mathtt{L}^{-\mathtt{T}}\mathbf{d}).

Figure 8. Example of Poisson’s summation formula.

Lemma 2.2.

Consider the distributions ρ(𝐫)=i=1Nbiδ(𝐫-𝐫i)ρ𝐫superscriptsubscripti1Nsubscriptbiδ𝐫subscript𝐫i\rho(\mathbf{r})=\sum_{i=1}^{N}b_{i}\,\delta(\mathbf{r}-\mathbf{r}_{i}), η(𝐫)=ρ(-𝐫)η𝐫ρ𝐫\eta(\mathbf{r})=\rho(-\mathbf{r}), and D(𝐫)=𝐧3δ(𝐫-𝙻𝐧)D𝐫subscript𝐧superscript3δ𝐫𝙻𝐧D(\mathbf{r})=\sum_{\mathbf{n}\in\mathbb{Z}^{3}}\delta(\mathbf{r}-\mathtt{L}% \mathbf{n}). The following identity holds:

𝐧3i=1Nj=1Nbibjψ(𝐫i-𝐫j+𝙻𝐧)=ψ(r)(ρηD)(𝐫)d𝐫,subscript𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j𝙻𝐧ψrnormal-⋆ρηD𝐫normal-d𝐫\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=𝐫rnorm𝐫r=\|\mathbf{r}\|.

Proof.

First,

ρη(𝐫)normal-⋆ρη𝐫\displaystyle\rho\star\eta\,(\mathbf{r}) =ρ(𝐬)η(𝐫-𝐬)d𝐬absentρ𝐬η𝐫𝐬normal-d𝐬\displaystyle=\int\rho(\mathbf{s})\,\eta(\mathbf{r}-\mathbf{s})\,\mathrm{d}% \mathbf{s}
=ρ(𝐬)ρ(𝐬-𝐫)d𝐬absentρ𝐬ρ𝐬𝐫normal-d𝐬\displaystyle=\int\rho(\mathbf{s})\,\rho(\mathbf{s}-\mathbf{r})\,\mathrm{d}% \mathbf{s}
=i=1Nj=1Nbibjδ(𝐬-𝐫i)δ(𝐫-𝐬-𝐫j)d𝐬absentsuperscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjδ𝐬subscript𝐫iδ𝐫𝐬subscript𝐫jnormal-d𝐬\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}
=i=1Nj=1Nbibjδ(𝐫-(𝐫i-𝐫j)).absentsuperscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjδ𝐫subscript𝐫isubscript𝐫j\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

ρηD=(ρη)D,normal-⋆ρηDnormal-⋆normal-⋆ρηD\rho\star\eta\star D=\left(\rho\star\eta\right)\star D,

so

ρηD(𝐫)normal-⋆ρηD𝐫\displaystyle\rho\star\eta\star D\,(\mathbf{r}) =ρη(𝐬)D(𝐬-𝐫)d𝐬absentnormal-⋆ρη𝐬D𝐬𝐫normal-d𝐬\displaystyle=\int\rho\star\eta\,(\mathbf{s})\,D(\mathbf{s}-\mathbf{r})\,% \mathrm{d}\mathbf{s}
=i=1Nj=1Nbibjδ(𝐬-(𝐫i-𝐫j))D(𝐬-𝐫)d𝐬absentsuperscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjδ𝐬subscript𝐫isubscript𝐫jD𝐬𝐫normal-d𝐬\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}
=i=1Nj=1NbibjD(𝐬-(𝐫i-𝐫j))absentsuperscriptsubscripti1Nsuperscriptsubscriptj1NsubscriptbisubscriptbjD𝐬subscript𝐫isubscript𝐫j\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}b_{i}b_{j}\,D(\mathbf{s}-(\mathbf{r}% _{i}-\mathbf{r}_{j}))
=𝐧3i=1Nj=1Nbibjδ(𝐬-(𝐫i-𝐫j+𝙻𝐧)).absentsubscript𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjδ𝐬subscript𝐫isubscript𝐫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,

ψ(r)(ρηD)(𝐫)d𝐫ψrnormal-⋆ρηD𝐫normal-d𝐫\displaystyle\int\psi(r)\,\left(\rho\star\eta\star D\right)(\mathbf{r})\,% \mathrm{d}\mathbf{r} =ψ(r)𝐧3i=1Nj=1Nbibjδ(𝐬-(𝐫i-𝐫j+𝙻𝐧))d𝐫absentψrsubscript𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjδ𝐬subscript𝐫isubscript𝐫j𝙻𝐧normal-d𝐫\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}
=𝐧3i=1Nj=1Nbibjψ(𝐫i-𝐫j+𝙻𝐧).absentsubscript𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j𝙻𝐧\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

U1=12ψ(r)(ρηD)(𝐫)d𝐫-12i=1Nbi2ψ(0),subscriptU112ψrnormal-⋆ρηD𝐫normal-d𝐫12superscriptsubscripti1Nsuperscriptsubscriptbi2ψ0U_{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 ρηDnormal-⋆ρηD\rho\star\eta\star D is real-valued, we can write

ψ(r)(ρηD)(𝐫)d𝐫=ψ^(h)(ρηD)^(𝐡)d𝐡,ψrnormal-⋆ρηD𝐫normal-d𝐫normal-^ψhnormal-^normal-⋆ρηD𝐡normal-d𝐡\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=𝐡hnorm𝐡h=\|\mathbf{h}\|. Consequently, equation (2.4) becomes

U1=12ψ^(h)(ρηD)^(𝐡)d𝐡-12i=1Nbi2ψ(0),subscriptU112normal-^ψhnormal-^normal-⋆ρηD𝐡normal-d𝐡12superscriptsubscripti1Nsuperscriptsubscriptbi2ψ0U_{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,

(ρηD)^(𝐡)=ρ^(𝐡)ρ^(-𝐡)1V𝐝3δ(𝐡-𝙻-𝚃𝐝)=|ρ^(𝐡)|21V𝐝3δ(𝐡-𝙻-𝚃𝐝).normal-^normal-⋆ρηD𝐡normal-^ρ𝐡normal-^ρ𝐡1Vsubscript𝐝superscript3δ𝐡superscript𝙻𝚃𝐝superscriptnormal-^ρ𝐡21Vsubscript𝐝superscript3δ𝐡superscript𝙻𝚃𝐝\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

ψ^(ρηD)^normal-^ψnormal-^normal-⋆ρηD\displaystyle\int\widehat{\psi}\,\widehat{\left(\rho\star\eta\star D\right)} =ψ^(h)(|ρ^(𝐡)|21V𝐝3δ(𝐡-𝙻-𝚃𝐝))d𝐡absentnormal-^ψhsuperscriptnormal-^ρ𝐡21Vsubscript𝐝superscript3δ𝐡superscript𝙻𝚃𝐝normal-d𝐡\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}
=1V𝐝3ψ^(𝙻-𝚃𝐝)|ρ^(𝙻-𝚃𝐝)|2absent1Vsubscript𝐝superscript3normal-^ψnormsuperscript𝙻𝚃𝐝superscriptnormal-^ρsuperscript𝙻𝚃𝐝2\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}
=1V𝐦ψ^(m)|ρ^(𝐦)|2absent1Vsubscript𝐦normal-^ψmsuperscriptnormal-^ρ𝐦2\displaystyle=\frac{1}{V}\sum_{\mathbf{m}}\widehat{\psi}(m)\,\big|\widehat{% \rho}(\mathbf{m})\big|^{2}
=1V𝐦0ψ^(m)|ρ^(𝐦)|2+1Vlimm0ψ^(m)|ρ^(0)|2absent1Vsubscript𝐦0normal-^ψmsuperscriptnormal-^ρ𝐦21Vsubscriptnormal-→m0normal-^ψmsuperscriptnormal-^ρ02\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}
=1V𝐦0ψ^(m)|ρ^(𝐦)|2+1Vlimm0ψ^(m)(i=1Nbi)2,absent1Vsubscript𝐦0normal-^ψmsuperscriptnormal-^ρ𝐦21Vsubscriptnormal-→m0normal-^ψmsuperscriptsuperscriptsubscripti1Nsubscriptbi2\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 𝐦=𝙻-𝚃𝐝𝐦superscript𝙻𝚃𝐝\mathbf{m}=\mathtt{L}^{-\mathtt{T}}\mathbf{d} and m=𝐦mnorm𝐦m=\|\mathbf{m}\|. So, to sum up, we have arrived at the identity:

U1subscriptU1\displaystyle U_{1} =12i=1Nijj=1Nbibjψ(𝐫i-𝐫j)+12𝐧0𝐧3i=1Nj=1Nbibjψ(𝐫i-𝐫j+𝙻𝐧)absent12superscriptsubscripti1Nsuperscriptsubscriptsuperscriptijj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j12subscriptsuperscript𝐧0𝐧superscript3superscriptsubscripti1Nsuperscriptsubscriptj1Nsubscriptbisubscriptbjψnormsubscript𝐫isubscript𝐫j𝙻𝐧\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}\|)
=12V𝐦0ψ^(m)|ρ^(𝐦)|2+12Vψ^(0)(i=1Nbi)2-12ψ(0)i=1Nbi2,absent12Vsubscript𝐦0normal-^ψmsuperscriptnormal-^ρ𝐦212Vnormal-^ψ0superscriptsuperscriptsubscripti1Nsubscriptbi212ψ0superscriptsubscripti1Nsuperscriptsubscriptbi2\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 𝒪(N3/2)𝒪superscriptN32\mathcal{O}(N^{3/2}) operations [dPS80] but it can be modified so that it involves only about 𝒪(NlogN)𝒪NN\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

𝐦0ψ^(m)|ρ^(𝐦)|2.subscript𝐦0normal-^ψmsuperscriptnormal-^ρ𝐦2\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

ρ^(𝐦)=j=1Nbje2πi𝐦𝚃𝐫j,normal-^ρ𝐦superscriptsubscriptj1Nsubscriptbjsuperscriptnormal-e2πnormal-isuperscript𝐦𝚃subscript𝐫j\widehat{\rho}(\mathbf{m})=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,% \mathbf{m}^{\mathtt{T}}\mathbf{r}_{j}},

where 𝐫j=𝙻𝐱jsubscript𝐫j𝙻subscript𝐱j\mathbf{r}_{j}=\mathtt{L}\,\mathbf{x}_{j} for some 𝐱j[0,1)3subscript𝐱jsuperscript013\mathbf{x}_{j}\in[0,1)^{3}. Moreover,

ρ^(𝐦)=j=1Nbje2πi𝐦𝚃𝐫j=j=1Nbje2πi𝐝𝚃𝙻-1𝙻𝐱j=j=1Nbje2πi𝐝𝚃𝐱j.normal-^ρ𝐦superscriptsubscriptj1Nsubscriptbjsuperscriptnormal-e2πnormal-isuperscript𝐦𝚃subscript𝐫jsuperscriptsubscriptj1Nsubscriptbjsuperscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙻1𝙻subscript𝐱jsuperscriptsubscriptj1Nsubscriptbjsuperscriptnormal-e2πnormal-isuperscript𝐝𝚃subscript𝐱j\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 ppp B-spline function is the ppp-fold convolution of the uniform density in the unit interval (see Figure 9).

Figure 9. Basic B-splines

By construction, each order ppp B-spline is of class Cp-1superscriptCp1C^{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 K1,K2,K3subscriptK1subscriptK2subscriptK3K_{1},K_{2},K_{3}\in\mathbb{N} and let

𝙺=[K1K2K3].𝙺delimited-[]subscriptK1absentabsentabsentsubscriptK2absentabsentabsentsubscriptK3\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

e2πi𝐦𝚃𝐫superscriptnormal-e2πnormal-isuperscript𝐦𝚃𝐫\displaystyle\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{m}^{\mathtt{T}}\mathbf{r}} =e2πi𝐝𝚃𝐱=e2πi𝐝𝚃𝙺-1𝐮3W(𝐮-)e2πi𝐝𝚃𝙺-1absentsuperscriptnormal-e2πnormal-isuperscript𝐝𝚃𝐱superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐮subscriptbold-ℓsuperscript3W𝐮bold-ℓsuperscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1bold-ℓ\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}}}
=𝐤𝐧3W(𝐮-(𝐤+𝙺𝐧))e2πi𝐝𝚃𝙺-1(𝐤+𝙺𝐧)absentsubscript𝐤subscript𝐧superscript3W𝐮𝐤𝙺𝐧superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺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}+\mathtt{K}\mathbf{n})}
=𝐤𝐧3W(𝐮-(𝐤+𝙺𝐧))e2πi𝐝𝚃𝙺-1𝐤e2πi𝐝𝚃𝐧=1absentsubscript𝐤subscript𝐧superscript3W𝐮𝐤𝙺𝐧superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐤subscriptnormal-⏟superscriptnormal-e2πnormal-isuperscript𝐝𝚃𝐧absent1\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}
=𝐤𝐧3W(𝐮-(𝐤+𝙺𝐧))e2πi𝐝𝚃𝙺-1𝐤,absentsubscript𝐤subscript𝐧superscript3W𝐮𝐤𝙺𝐧superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺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

=𝐤+𝙺𝐧withki{0,1,,Ki-1}and𝐧3.formulae-sequencebold-ℓ𝐤𝙺𝐧withformulae-sequencesubscriptki01normal-…subscriptKi1and𝐧superscript3\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𝐤subscriptQ𝐤Q_{\mathbf{k}} as

Q^𝐝=DFT[Q𝐤]𝐝=k1=0K1-1k2=0K2-1k3=0K3-1Q𝐤e-2πi𝐝𝚃𝙺-1𝐤.subscriptnormal-^Q𝐝DFTsubscriptdelimited-[]subscriptQ𝐤𝐝superscriptsubscriptsubscriptk10subscriptK11superscriptsubscriptsubscriptk20subscriptK21superscriptsubscriptsubscriptk30subscriptK31subscriptQ𝐤superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐤\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

IDFT[Q^𝐝]𝐤=d1=0K1-1d2=0K2-1d3=0K3-1Q^𝐝e2πi𝐝𝚃𝙺-1𝐤.IDFTsubscriptdelimited-[]subscriptnormal-^Q𝐝𝐤superscriptsubscriptsubscriptd10subscriptK11superscriptsubscriptsubscriptd20subscriptK21superscriptsubscriptsubscriptd30subscriptK31subscriptnormal-^Q𝐝superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐤\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

ρ^(𝐦)normal-^ρ𝐦\displaystyle\widehat{\rho}(\mathbf{m}) =j=1Nbje2πi𝐝𝚃𝐱jabsentsuperscriptsubscriptj1Nsubscriptbjsuperscriptnormal-e2πnormal-isuperscript𝐝𝚃subscript𝐱j\displaystyle=\sum_{j=1}^{N}b_{j}\,\mathrm{e}^{2\pi\mathrm{i}\,\mathbf{d}^{% \mathtt{T}}\mathbf{x}_{j}}
𝐤j=1Nbj𝐧3W(𝐮-(𝐤+𝙺𝐧))=Q𝐤e2πi𝐝𝚃𝙺-1𝐤absentsubscript𝐤subscriptnormal-⏟superscriptsubscriptj1Nsubscriptbjsubscript𝐧superscript3W𝐮𝐤𝙺𝐧absentsubscriptQ𝐤superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐤\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}}
=𝐤Q𝐤e2πi𝐝𝚃𝙺-1𝐤=Q^𝐝.absentsubscript𝐤subscriptQ𝐤superscriptnormal-e2πnormal-isuperscript𝐝𝚃superscript𝙺1𝐤subscriptnormal-^Q𝐝\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:

ρQ𝐼𝐹𝑇𝐼𝐷𝐹𝑇ρ^Q^ρsuperscriptnormal-→Qnormal-↓ IFT absentmissing-subexpressionnormal-↓absent IDFT missing-subexpressionmissing-subexpressionnormal-^ρsuperscriptnormal-→normal-^Q\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:

𝐦0ψ^(m)|ρ^(𝐦)|2=𝐦0ψ^(m)ρ^(-𝐦)ρ^(𝐦)𝐝0ψ^(𝙻-𝚃𝐝)Q^-𝐝Q^𝐝subscript𝐦0normal-^ψmsuperscriptnormal-^ρ𝐦2subscript𝐦0normal-^ψmnormal-^ρ𝐦normal-^ρ𝐦subscript𝐝0normal-^ψnormsuperscript𝙻𝚃𝐝subscriptnormal-^Q𝐝subscriptnormal-^Q𝐝\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

A𝐝={φ^(𝙻𝚃𝐡),if 𝐡0𝐡0\mathbf{h}\neq 0,0,if 𝐡=0𝐡0\mathbf{h}=0.subscriptA𝐝casesnormal-^φnormsuperscript𝙻𝚃𝐡if 𝐡≠0𝐡0\mathbf{h}\neq 00if 𝐡=0𝐡0\mathbf{h}=0A_{\mathbf{d}}=\begin{cases}\widehat{\varphi}(\|\mathtt{L}^{\mathtt{T}}\,% \mathbf{h}\|),&\text{if $\mathbf{h}\neq 0$},\\ 0,&\text{if $\mathbf{h}=0$}.\end{cases}

with

hi(𝐝)={di,if 0di<12Ki0subscriptdi12subscriptKi0\leq d_{i}<\frac{1}{2}K_{i},Ki-di,if 12Kidi<Ki12subscriptKisubscriptdisubscriptKi\frac{1}{2}K_{i}\leq d_{i}<K_{i}.subscripthi𝐝casessubscriptdiif 0≤di<12⁢Ki0subscriptdi12subscriptKi0\leq d_{i}<\frac{1}{2}K_{i}subscriptKisubscriptdiif 12⁢Ki≤di<Ki12subscriptKisubscriptdisubscriptKi\frac{1}{2}K_{i}\leq d_{i}<K_{i}h_{i}(\mathbf{d})=\begin{cases}d_{i},&\text{if $0\leq d_{i}<\frac{1}{2}K_{i}$}% ,\\ K_{i}-d_{i},&\text{if $\frac{1}{2}K_{i}\leq d_{i}<K_{i}$}.\end{cases}

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

d1=-K12K12-1d2=-K22K22-1d3=-K32K32-1ψ^(𝙻-𝚃𝐝)Q^-𝐝Q^𝐝=d1=0K1-1d2=0K2-1d3=0K3-1A𝐝Q^-𝐝Q^𝐝.superscriptsubscriptsubscriptd1subscriptK12subscriptK121superscriptsubscriptsubscriptd2subscriptK22subscriptK221superscriptsubscriptsubscriptd3subscriptK32subscriptK321normal-^ψnormsuperscript𝙻𝚃𝐝subscriptnormal-^Q𝐝subscriptnormal-^Q𝐝superscriptsubscriptsubscriptd10subscriptK11superscriptsubscriptsubscriptd20subscriptK21superscriptsubscriptsubscriptd30subscriptK31subscriptA𝐝subscriptnormal-^Q𝐝subscriptnormal-^Q𝐝\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]),

f^dg=g^df,normal-^fnormal-dgnormal-^gnormal-df\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)

ψ^η^ρ^=ψ^η^^ρ=(C^ρ)ρnormal-^ψnormal-^ηnormal-^ρnormal-^normal-^ψnormal-^ηρnormal-⋆normal-^Cρρ\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 fg^=f^g^normal-^fgnormal-⋆normal-^fnormal-^g\widehat{fg}=\widehat{f}\star\widehat{g} and that, for f=f(x)ffxf=f(x), f^^(x)=f(-x)normal-^normal-^fxfx\widehat{\widehat{f}}(x)=f(-x) (again, see [WOL03]).

Therefore, the algorithm boils down to:

Q𝐤subscriptQ𝐤\displaystyle Q_{\mathbf{k}} Q^-𝐝=DFT(Q)dnormal-→absentsubscriptnormal-^Q𝐝DFTsubscriptQd\displaystyle\rightarrow\widehat{Q}_{-\mathbf{d}}=\text{DFT}(Q)_{d}
A𝐝DFT(Q)𝐝normal-→absentnormal-⋅subscriptA𝐝DFTsubscriptQ𝐝\displaystyle\rightarrow A_{\mathbf{d}}\cdot\text{DFT}(Q)_{\mathbf{d}}
IDFT(ADFT(Q))𝐤=(IDFT(A)Q)𝐤normal-→absentIDFTsubscriptnormal-⋅ADFTQ𝐤subscriptnormal-⋆IDFTAQ𝐤\displaystyle\rightarrow\text{IDFT}(A\cdot\text{DFT}(Q))_{\mathbf{k}}=(\text{% IDFT}(A)\star Q)_{\mathbf{k}}
(IDFT(A)Q)𝐤Q𝐤,normal-→absentnormal-⋅subscriptnormal-⋆IDFTAQ𝐤subscriptQ𝐤\displaystyle\rightarrow(\text{IDFT}(A)\star Q)_{\mathbf{k}}\cdot Q_{\mathbf{k% }},

Appendix A Assorted results in Harmonic analysis

Theorem A.1 (The convolution is commutative).

Theorem A.2 (The Fourier transform of a convolution is the product of the Fourier transforms).

Theorem A.3 ([WOL03]).

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

μ^dν=ν^dμnormal-^μnormal-dνnormal-^νnormal-dμ\int\hat{\mu}\,\mathrm{d}\nu=\int\hat{\nu}\,\mathrm{d}\mu

Theorem A.4 (Plancherel’s theorem).

References

  • [DYP93] T. Darden, D. York and L. Pedersen (1993) Particle mesh Ewald: An Nlog(N)NNN\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-nsuperscriptRn{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.