Elsevier

Journal of Computational Physics

Volume 347, 15 October 2017, Pages 341-366
Journal of Computational Physics

The Spectral Ewald method for singly periodic domains

https://doi.org/10.1016/j.jcp.2017.07.001Get rights and content

Abstract

We present a fast and spectrally accurate method for efficient computation of the three dimensional Coulomb potential with periodicity in one direction. The algorithm is FFT-based and uses the so-called Ewald decomposition, which is naturally most efficient for the triply periodic case. In this paper, we show how to extend the triply periodic Spectral Ewald method to the singly periodic case, such that the cost of computing the singly periodic potential is only marginally larger than the cost of computing the potential for the corresponding triply periodic system.

In the Fourier space contribution of the Ewald decomposition, a Fourier series is obtained in the periodic direction with a Fourier integral over the non-periodic directions for each discrete wave number. We show that upsampling to resolve the integral is only needed for modes with small wave numbers. For the zero wave number, this Fourier integral has a singularity. For this mode, we effectively need to solve a free-space Poisson equation in two dimensions. A very recent idea by Vico et al. makes it possible to use FFTs to solve this problem, allowing us to unify the treatment of all modes. An adaptive 3D FFT can be established to apply different upsampling rates locally. The computational cost for other parts of the algorithm is essentially unchanged as compared to the triply periodic case, in total yielding only a small increase in both computational cost and memory usage for this singly periodic case.

Introduction

In molecular dynamics simulations, a crucial and time consuming task is to compute the long-range interactions or particularly the electrostatic potential (or force) between charged particles. For systems that are subject to periodic boundary conditions, Ewald summation is a technique to evaluate these interactions. For N particles with charges qn at positions xnXi=13[0,Li), n=1,,N, the electrostatic potential or briefly the potential evaluated at a target point xm is written asφ(xm)=pPDn=1Nqn|xmxn+p|, where PD with D=1,2,3, can be modified such that it expresses the periodicity. The prime denotes that the term with n=m is omitted from the sum for p=0. To impose the periodicity in three dimensions (3d-periodicity), the simulation box is replicated in three directions and we define P3={(α1L1,α2L2,α3L3):αiZ}. In this case, in light of the neutrality condition, i.e., n=1Nqn=0, the sum in (1.1) is only conditionally convergent and therefore the order of summation has to be exactly defined. In the classical Ewald sum proposed by Ewald [14], (1.1) is decomposed into a fast decaying part and a smooth part which is computed in Fourier space. The result is that of a spherical summation order, and the sum can be written asφ3P(xm)=φmR+φmF+φmself=pP3n=1Nqnerfc(ξ|xmxn+p|)|xmxn+p|+4πVk0ek2/4ξ2k2n=1Nqneik(xmxn)2ξπqm, where m=1,,N, k=2π(n1L1,n2L2,n3L3), with niZ and k=|k|. The superscripts R, F, and self denote the real-space, Fourier-space (here referred to as k-space), and self correction term, respectively. The self correction term is added in order to eliminate the self interaction contribution of the charges included due to the decomposition. Moreover, applying the spherical order of summation, the k=0 term (dipole term) of the 3d periodic Ewald sum, depends on the dielectric constant of the surrounding medium. If the medium has an infinite dielectric constant, the dipole term vanishes, cf. [15]. Assuming so, the k=0 term is excluded from the sum.

In (1.2), ξ>0 is the decomposition parameter (Ewald parameter) which controls the decay of the terms in the real space and k-space sums, but does not alter the total result. Under the assumption of a uniform distribution of charges, a proper choice of ξ can decrease the computational complexity of computing the potential (1.2) at xm, m=1,,N, from O(N2) to O(N3/2) albeit with a very large constant. This however can be reduced to O(Nlog(N)) (also with a much smaller constant) using fast methods which take advantage of the Fast Fourier transform (FFT). Inspired by the Particle–Particle–Particle Mesh Ewald (P3M) method by Hockney and Eastwood [17], different fast methods have been proposed, including the Smooth Particle mesh Ewald (SPME) method [13], and a spectrally accurate Ewald method for triply periodic (SE3P) [19] and doubly periodic (SE2P) [20] systems. Using the idea in [30] we have developed another spectrally accurate method for the fast evaluation of sums involving free-space Green's functions [1].

In the present work, we complete the framework of the spectral Ewald methods by extending the algorithm to systems with one periodic direction in three dimensions (1d- or singly periodic). This method can be used, e.g., for simulation of nanopores [6] and nanotubes [9]. In one of the first attempts, Lekner summation was used to compute the long-range interactions in 1d-periodic system of particles [7], [3], [8]. The first derivation of the 1d-periodic Ewald summation Ewald1P was given by Porto [23] with an integral representation of the k-space sum. The integral can however be evaluated to obtain the closed form of the formula, see [27, Appendix D] and references therein. Recently, Nestler et al. [22] developed a fast algorithm which employs non-equispaced FFTs (NFFT). To the best of our knowledge, this is the only method with O(Nlog(N)) complexity for 1d-periodic problems. The approach that we have taken differs significantly from theirs, as will be commented on in Section 3.

FFT based methods such as methods in the PME (particle mesh Ewald) family are most efficient for the triply periodic case. In this case, FFTs can be used in all directions without any oversampling. As soon as there is a non-periodic direction, the grid has to be extended in that direction. In the doubly periodic case, Arnold et al. [4], [10] devised a method where the problem is extended to full periodicity, with a larger length in the non-periodic direction, and where a correction term is applied to improve on the result. Here, the increased length in the non-periodic direction simply means a zero-padding of the FFT, increasing in the number of grid points in that direction. The SE2P method by Lindbo and Tornberg [20] takes a different approach, which needs a “mixed” transform; a discrete Fourier transform in the periodic variables and an approximation to the continuous Fourier integral transform in the free dimension. Also in this case the grid in the free dimension must be oversampled for an accurate approximation.

Extending the Spectral Ewald method to the singly periodic case, there were two main challenges to overcome. Firstly, the oversampling. An oversampling by a factor of four to six makes the FFTs four to six times more expensive to compute when the oversampling is applied in one dimension. With two free dimensions, this would increase the cost by a factor of 16 to 64, which is clearly not desirable. Oversampling needs to be done to resolve the Fourier integrals. By analysis of a similar one-dimensional integral we can understand how the error behaves and recognize that only for small discrete wave numbers (a small number of periodic modes) do the FFT grids need to be upsampled. Based on this, we have developed what we call adaptive FFTs and IFFTs (denoted by AFT and AIFT in this paper) that only upsample for a select number of discrete modes in the periodic direction. The ratio of the run time for the AFT and the FFT without oversampling decreases with grid size, as a smaller fraction of modes must be oversampled, and a typical increase in cost can be a factor of 2–3 instead.

In the derivation of the singly periodic Ewald sums, there is a term that includes the contribution from the zero wave number in the periodic direction, i.e., that depends only on the variables in the free directions. The direct evaluation of this sum at all target points would however incur an O(N2) computational cost, and the second challenge was to significantly reduce this. One interpretation of this sum is that it is the solution to the Poisson equation in R2, with the right hand side convolved with specifically scaled Gaussians centered at each of the charge locations and projected onto the plane z=0, with z the periodic direction. A very recent idea for how to solve free space problems by the means of FFTs [30] can therefore be used. Hence, we are able to include this zero wave number contribution into the treatment of the full k-space term. This extra sum then only entails a special scaling in Fourier space for the modes with a zero wave number in the periodic direction. This is done at a negligible extra cost.

With these two main advances, we have developed a fast and spectrally accurate FFT-based method for the evaluation of the k-space sum in the Ewald summation formula for the singly periodic case. We will denote this method the SE1P method.

The outline of this paper is as follows; In section 2, we present the Ewald1P formulas and provide the truncation error estimates. We introduce our fast method for computing the electrostatic force and potential in section 3. Section 4 is devoted to a discussion regarding approximation errors including the errors introduced by discretization of the Fourier integrals, and the related issue of parameter selection including oversampling rates in the adaptive FFT. The following section is dedicated to the implementation details of the algorithm. Finally, we supply numerical results including comparisons with the 3d-periodic counterpart in section 6 and wrap up with conclusions in section 7.

Section snippets

The Ewald summation formula

In this section, we present the three dimensional Ewald summation formula under 1d-periodic boundary conditions. The formulas to compute (1.1) can be derived using Fourier integrals. The reader may consult [27] for details and alternatives of the derivation. We note here that unlike the 3d-periodic case, equation (1.1) is shown to be absolutely convergent in the 1d-periodic case under the assumption of charge neutrality and is therefore independent of the summation order [3].

Henceforth we

Introducing a fast method

According to the discussion in the previous section, contributions to the real space sum will be ignored if the distance between the source location and the target evaluation point is larger than a cut-off radius rc. Typically, a linked cell list [2], [17] or a Verlet list algorithm [2], [29] can be used to efficiently obtain a list of nearest neighbors. The real space sum in (2.2) includes a summation over the one periodic dimension, which means that contribution from periodic images of the

Approximation errors and parameter selection

In addition to the truncation error due to the finite representation of the k-space sum (section 2.2), approximation errors are also involved in our fast method. These errors are committed due to (a) the truncation of Gaussians, (b) applying the quadrature rule to evaluate (3.24) and (c) approximating Fourier integrals. The approximation errors due to (a) and (b) are considered in section 4.1 and (c) is discussed in section 4.2.

AFT/AIFT

We now present an algorithm which accelerates the computation of mixed Fourier transforms introduced in [20]. As has been stated before, our system of interest is periodic in z and free in the x and y directions. Consider again the local pad set I (3.2), and the associated set J (3.3), local pad size nl and oversampling factors sl and s0 introduced in section 3.2. Moreover, let M be the number of grid points and let k32πL3{M2,,M21} and κ1,κ22πsL˜{M˜2,,M˜21}, with s defined in (3.4). The

Numerical results

In this section we present numerical results of computing the electrostatic potential and force with 1d-periodicity using the SE1P method. All the simulations are done on one core on a machine with Intel Core i7-3770 CPU which runs on 3.40 GHz with 8 GB of memory. The FFT/IFFT and scaling steps are done in MATLAB and the gridding and interpolation steps are written in C and are dynamically linked and called through the MATLAB MEX interface. The subroutines are written in C and are built with

Summary and conclusions

We develop a fast and accurate algorithm to compute the electrostatic potential, force and energy for three dimensional systems of charged particles under singly periodic boundary conditions. The method is based on the Ewald summation formula and follows the general framework of other Particle mesh Ewald (PME) methods with an FFT treatment of the Fourier sum. Specifically, this work is an extension of the Spectral Ewald method that has been developed for triply periodic (SE3P, [19]) and doubly

Acknowledgements

This work has been supported by the Swedish Research Council under grant no. 2011-3178 and by the Swedish e-Science Research Center. The authors gratefully acknowledge this support.

References (31)

  • A. Arnold et al.

    Comparison of scalable fast methods for long-range interactions

    Phys. Rev. E

    (2013)
  • I.C. Bourg et al.

    Molecular dynamics simulations of water structure and diffusion in silica nanopores

    J. Phys. Chem. C

    (2012)
  • A. Bródka et al.

    Electrostatic interactions in computer simulations of a three-dimensional system periodic in two directions: Ewald-type summation

    J. Chem. Phys.

    (2002)
  • A. Bródka et al.

    Three-dimensional Ewald method with correction term for a system periodic in one direction

    J. Chem. Phys.

    (2004)
  • J. De Joannis et al.

    Electrostatics in periodic slab geometries. II

    J. Chem. Phys.

    (2002)
  • Cited by (0)

    View full text