A fast spectral method for electrostatics in doubly-periodic slit channels
Ondrej Maxian, Raul P. Peláez, Leslie Greengard, Aleksandar Donev
AA fast spectral method for electrostatics in doubly-periodic slit channels
Ondrej Maxian, Raul P. Pel´aez, Leslie Greengard,
1, 3 and Aleksandar Donev Courant Institute, NYU, New York, NY 10012 Department of Theoretical Condensed Matter Physics,Universidad Aut´onoma de Madrid, 28049, Madrid, Spain Center for Computational Mathematics,Flatiron Institute, New York, NY 10010
We develop a fast method for computing the electrostatic energy and forces for acollection of charges in doubly-periodic slabs with jumps in the dielectric permittiv-ity at the slab boundaries. Our method achieves spectral accuracy by using Ewaldsplitting to replace the original Poisson equation for nearly-singular sources with asmooth far-field Poisson equation, combined with a localized near-field correction.Unlike existing spectral Ewald methods, which make use of the Fourier transformin the aperiodic direction, we recast the problem as a two-point boundary valueproblem in the aperiodic direction for each transverse Fourier mode, for which exactanalytic boundary conditions are available. We solve each of these boundary valueproblems using a fast, well-conditioned Chebyshev method. In the presence of dielec-tric jumps, combining Ewald splitting with the classical method of images results insmoothed charge distributions which overlap the dielectric boundaries themselves.We show how to preserve high order accuracy in this case through the use of aharmonic correction which involves solving a simple Laplace equation with smoothboundary data. We implement our method on Graphical Processing Units, andcombine our doubly-periodic Poisson solver with Brownian Dynamics to study theequilibrium structure of double layers in binary electrolytes confined by dielectricboundaries. Consistent with prior studies, we find strong charge depletion near theinterfaces due to repulsive interactions with image charges, which points to the needfor incorporating polarization effects in understanding confined electrolytes, boththeoretically and computationally. a r X i v : . [ m a t h . NA ] J a n I. INTRODUCTION
The evaluation of electrostatic interactions in a collection of charges is a classical prob-lem in computational physics, with applications to the study of electrolyte solutions, macro-molecules, ion channels, and other systems. In Molecular Dynamics (MD) and BrownianDynamics (BD) methods, the forces between the charges need to be computed at least onceper time step, while in Monte Carlo (MC) methods the total electrostatic energy of the entirecollection of charges is required. Hybrid BD-MC schemes require computing both forces andenergy. Each of the many linear-scaling methods available to compute these quantities (forvarying geometries and boundary conditions) falls into one of two categories: fast multipolemethods [1–3], and variants of particle-(particle-particle)-mesh (P3M) methods, includingpre-corrected FFT and spectral Ewald (SE) methods [4–9]. Our focus in this work is on thelatter type of method, which tends to be more efficient than the former for homogeneouscharge distributions.In this paper we develop a linear-scaling spectral Ewald method to compute electro-static forces and energy in doubly-periodic or slab geometries with dielectric jumps at theboundaries of the slab, which is applicable to the study of (Debye) double layers in confinedelectrolyte solutions. Specifically, we will consider charges in a slab 0 < z < H immersedin a uniform dielectric medium in a domain that is periodic in the xy plane but potentiallyhas jumps in the dielectric permittivity at z = 0 and z = H . We will consider not pointcharges, but rather Gaussian charges of width g w >
0. This approach avoids divergentself-interactions and unnecessarily stiff electrostatic interactions at short distances and isconsistent with our focus on electrolyte solutions, for which the charged particles are sol-vated ions that are not actually point charges to begin with. Nevertheless, our methodapplies to arbitrarily small g w (to within roundoff errors), allowing us to approach the limitof point charges if desired.There are already a large number of methods in the literature for doubly-periodic elec-trostatics [3–5, 10–14]. We will not attempt to review and compare all of them, but focusinstead on highlighting the improvements of our approach over recent SE methods [5, 14],to which our approach is most closely related. Our method combines a number of ideasfrom the existing literature with some new ideas and new numerical methods, which resultsin an improved P3M method in the slab geometry. Like many previous approaches, we usethe method of images to tackle the presence of dielectric jumps [3, 10, 13] and convert theproblem to a doubly-periodic problem in a medium with uniform dielectric permittivity. Atthe same time, however, our method avoids the need to consider the full (infinite) imagesystem by imposing exact boundary conditions at z = 0 and z = H , and using a grid-based Poisson/Laplace solver to account for the distant images. We make key use of Ewaldsplitting, so that some aspects of our method are very similar to existing SE approaches.Our use of Ewald splitting is different from that in [5, 14]. The traditional view is toconsider Ewald splitting as an analytical technique that separates the electrostatic sumsfor point-like charges into “real-space” and “Fourier-space” contributions. The real-spaceor near-field part is easy to handle by direct summation over pairs of nearby particles,while the Fourier-space or far-field part is handled using Fast Fourier Transform (FFT)based methods. These are simplest to understand in triply-periodic domains [6], and thecorresponding formulas are easy to derive. This standard view has been extended to othergeometries including slabs, as reviewed in [7, 15].One drawback of the traditional view of Ewald splitting is that, by focusing on the Fourierdomain alone, non-physical sampling requirements are sometimes imposed on the method.The paper [14], which was perhaps the first to introduce and carefully analyze a spectrallyaccurate approach for doubly-periodic geometries, requires large oversampling factors thatdepend on both the tolerance and the aspect ratio of the domain. While a later paper[5] was able to overcome this, we believe the applicability of Ewald-type methods to moregeneral geometries is easier to understand in the framework we present here, which makessimultaneous use of both Fourier and more standard PDE-based ideas. This framework iskey to understanding the novel parts of our algorithm.Thus, we begin by considering Ewald splitting as an improvement of P3M approaches,where the first step is to create a smooth source distribution by convolution with a Gaussian,the second step involves solving the Poisson equation with this smooth right-hand side, andthe third step is to correct for the (localized) errors introduced by the initial smoothing. Onecan use any sufficiently accurate solver for the smoothed problem, not necessarily one basedon Fourier analysis. The third step is handled (as in any Ewald-type method) analytically,using pairwise summation over some collection of near neighbors.Having split the problem into a PDE with smooth data (accurately representing thefar-field interactions) and a near-field correction, it remains only to correctly specify appro-priate boundary conditions in the unbounded z direction and to construct an appropriatefast solver. It is straightforward to derive a Robin-type condition on the slab boundary foreach mode, through an analysis of the Dirichlet-to-Neumann (DtN) map. After Fouriertransformation with respect to the periodic directions, we use a real-space Chebyshev spec-tral solver in the now finite z direction. Unlike pure Fourier-based schemes, such as [14],all that is required is that the source distribution be resolved by the Chebyshev grid. Nooversampling is required and the aspect ratio of the domain plays no role. Unlike [5], wemaintain spectral accuracy rather than switching over to an algebraic convergence rate.That said, since Chebyshev methods are easily implemented efficiently and robustly usinga Fourier transform on a double-sized grid, we still make use of three-dimensional FFTs asthe key component to achieving O ( N log N ) scaling where N is the number of unknowns.The fixed oversampling factor of 2 is typical for any FFT-based aperiodic convolution.For a single dielectric interface, our Ewald splitting approach can be applied after usingthe classical image construction to handle the dielectric jump in the z direction. In particular,a set of images can be constructed and Ewald splitting can be used to smear the charges,thereby allowing for a coarse grid in the doubly-periodic Poisson solver. This approachruns into difficulties, however, for slabs with two (or more) dielectric jumps. Our methodavoids the inherent problem of having infinitely many images for a slab geometry by onlyincluding the images of charges that are sufficiently close to the dielectric boundary for thecorresponding Ewald cloud to overlap the boundary. We account for the rest of the infiniteimage system through a harmonic correction potential that can be computed analyticallyfrom the mismatch in the boundary conditions at the slab boundaries. This simple butpowerful idea appears have been overlooked in the field, and is easily combined with otherdoubly-periodic Poisson solvers, such as the nonuniform FFT-based solver proposed in [5].An outline of the paper follows. In Section II, we present the mathematical formulation ofthe problem. In Section III A, we develop a novel continuum approach for smooth doubly-periodic electrostatics problems based on the DtN map. We then present our variant ofEwald splitting in Section III B that maps a nonsmooth problem into a smooth one. InSection IV, we present our main contribution: a continuum approach to doubly-periodicelectrostatics for slabs that combines a restricted image construction, Ewald splitting, the The use of the DtN map is not new; see, for example, Eqs. (17,18) in [4] or the FMM-based algorithm in[3], but our use of it within the SE framework appears not to have been explored.
DtN-based doubly-periodic solver, and a correction approach for coarse-graining imagesthat are sufficiently far away from the slab. In Section V, we present a discretization ofthe continuum formulation using FFTs in the xy plane and Chebyshev polynomials in the z direction, which we implement in a public-domain code running on Graphical ProcessingUnits (GPUs). In Section VI, we validate the accuracy of our method by comparing toreference analytical and numerical results, and in Section VII we combine our electrostaticsolver with Brownian Dynamics to study binary electrolyte solutions in slit channels witheither uncharged (Section VII B) or charged walls (Section VII C). By comparing our resultsto reference Monte Carlo results from the literature and analytical solutions, we validateour method and establish the importance of polarization effects that come from jumps indielectric permittivity. We conclude with a summary and a discussion of future directionsin Section VIII. II. PROBLEM STATEMENT
We consider solving an electrostatics problem for a collection of N Gaussian charges withstrengths q k and positions z k . The domain geometry is that of a slit channel: periodic inthe x and y directions on [0 , L x ] and [0 , L y ], respectively, and unbounded in z . We willalso assume that the charges are contained within a finite region z ∈ (0 , H ), and that thereare fixed surface charge densities σ b/t ( x, y ) on the bottom and top boundaries of this region( z = 0 and z = H ). These assumptions give the electrostatic equation for the potential φ ( x = ( x, y, z )), − ∇ · ( (cid:101) (cid:15) ∇ φ ) = ρ, where (1) ρ ( x ) = ρ ( x, y, z ) = N (cid:88) k =1 q k g ( (cid:107) x − z k (cid:107) ) + σ b ( x, y ) δ ( z ) + σ t ( x, y ) δ ( z − H ) , (2)on ( x, y ) ∈ [0 , L x ] × [0 , L y ] with z unbounded. We assume that each charge has a Gaussiancharge density g ( r ) = 1(2 πg w ) / exp (cid:18) − r g w (cid:19) , (3)with standard deviation g w related to the physical size of the charges, with g w → φ is only unique up to a constant, and sowe set φ ( ) = 0 . (4)Our goal is to solve (1) for a slab with piecewise constant dielectric permittivity in the z direction, (cid:101) (cid:15) ( x ) = (cid:101) (cid:15) ( z ) = (cid:15) b z < (cid:15) < z < H(cid:15) t z > H (5)Substituting this into the electrostatic equation (1), we obtain a Poisson equation for z ∈ (0 , H ), (cid:15) ∆ φ ( x ) = − N (cid:88) k =1 q k g ( (cid:107) x − z k (cid:107) ) := − f ( x ) (6)together with boundary conditions on the potential and electric displacement at z = 0 and z = H , φ ( x, y, z → + ) = φ ( x, y, z → − ) (7) φ ( x, y, z → H − ) = φ ( x, y, z → H + ) (8) (cid:15) ∂φ∂z ( x, y, z → + ) − (cid:15) b ∂φ∂z ( x, y, z → − ) + σ b ( x, y ) = 0 , (9) (cid:15) ∂φ∂z ( x, y, z → H − ) − (cid:15) t ∂φ∂z ( x, y, z → H + ) − σ t ( x, y ) = 0 . (10)We assume here that the domain is overall electroneutral (including the wall-bound chargedensities), N (cid:88) k =1 q k + (cid:90) L x (cid:90) L y ( σ b ( x, y ) + σ t ( x, y )) dx dy = 0 , (11)and so the electric field ∇ φ must decay to zero as z → ±∞ . In (11), we have assumed thateach Gaussian charge density (3) is fully contained inside the slab, so that it integrates tounity on [0 , H ]. Of course, this is not possible exactly since a Gaussian is not compactlysupported. To address this, we truncate the Gaussian at a finite distance n σ g w ≥ g w , sothat the three-dimensional integral of each charge density is at least 99 .
9% of the charge q . For the Gaussian to be fully contained inside the slab, we assume that the truncatedGaussian envelopes do not overlap the dielectric boundaries, z k ∈ [ n σ g w , H − n σ g w ].For MCMC simulations, we need to compute the electrostatic energy U = 12 (cid:90) L x (cid:90) L y (cid:90) ∞−∞ φ ( x, y, z ) ρ ( x, y, z ) dz dy dx. (12)Substituting the expression (2) for the charge density ρ , we obtain the energy U = 12 N (cid:88) k =1 q k φ ( z k ) + 12 (cid:90) L y (cid:90) L x ( σ ( x, y, φ ( x, y,
0) + σ ( x, y, H ) φ ( x, y, H )) dx dy, (13)where φ ( z ) = (cid:90) x φ ( x ) g ( (cid:107) x − z (cid:107) ) d x (14)is the convolution of the pointwise potential φ ( x ) with a Gaussian of width g w . This meansthat the energy (13) is the sum of the average potential at the center of each charge plusthe energy due to the surface charge densities σ b/t .The electrostatic forces F k = − ∂U/∂ z k can now be determined in a straightforward wayfrom the energy (13). The pointwise electric field can be determined from the pointwisepotential by E = − ∇ φ, (15)and the average electric field is the convolution E ( z ) = (cid:90) x E ( x ) g ( (cid:107) x − z (cid:107) ) d x . (16)The force on each charge is given by F k = − ∂U∂ z k = q k E ( z k ) . (17)Our goal will be to compute U ( z ) and F ( z ) to high accuracy in (log) linear time in thenumber of charges N . III. SOLVER COMPONENTS
In this section, we build on prior work to introduce solvers for two simplified problems:the Poisson equation (6) for smooth charge density f and jump BCs (7) − (10), and a freespace solver for point charges. In Section IV, we combine these two pieces to yield a solverfor (6) − (10) that remains efficient as g w → smooth doubly periodic problems, wherea grid-based solver can efficiently resolve the charge density f ( x ), as would be the case forlarge g w . Our method in this case is to use superposition to split the problem into twosub-problems: a Poisson equation over free space with uniform permittivity, and a harmoniccorrection to account for the jump conditions. This solver is inefficient for (near) pointcharges (small g w ) because the number of grid points necessary to resolve f ( x ) becomes toolarge. We address this in Section III B via Ewald splitting. A. Smooth doubly periodic problems
Let us suppose first that the distance between the charges is comparable to their width g w . In this case, a grid-based method is an efficient way to solve (6) with the BCs (7) − (10).While this can be done with a single solve, it will aid us in Section IV to split the potentialinto two pieces, φ = φ ∗ + φ ( c ) , (18)and solve for each piece separately. The first of these, φ ∗ , is smooth and found by solvingthe Poisson equation (6), but with free space BCs and a uniform permittivity (cid:15) (see III A 1).The second piece, φ ( c ) , is found by solving a Laplace equation for the correction potentialthat satisfies the BCs (7) − (10) (see III A 2). This part of the solution is necessarily onlypiecewise smooth, since it is built to satisfy jump BCs. Because each of the pieces accountsfor either the Gaussian or wall-bound charge densities (but not both), neither of them isnecessarily electroneutral, and we therefore conclude in Section III A 3 by constructing acombined solution φ that gives a vanishing electric field as z → ±∞ .
1. First problem: Poisson solve in free space
For our first problem, we let φ ∗ be defined everywhere as the solution of a Poisson equationwith uniform permittivity (cid:15) and no wall-bound charge densities, (cid:15) ∆ φ ∗ ( x ) = − M (cid:88) k =1 q k g ( (cid:107) x − z k (cid:107) ) . (19)We take free space boundary conditions in z , so that ∂φ/∂z ( z → ±∞ ) →
0, and periodicboundary conditions in x and y on [0 , L x ] × [0 , L y ]. For the moment we have no recourseto a grid-based solver, since the z domain is unbounded. But since f is a sum of chargedensities concentrated on (0 , H ), the Poisson equation (19) is actually a Laplace equationin the exterior, (cid:15) ∆ φ ∗ = 0 on z < z > H. (20)This Laplace equation can be solved analytically by taking a Fourier transform in x and y ,which we denote with a hat. Using wave numbers k (cid:107) = ( k x , k y ) = (2 πn/L x , πm/L y ), with m and n integers and k (cid:107) = (cid:13)(cid:13) k (cid:107) (cid:13)(cid:13) , we obtain ∂ (cid:98) φ ∗ ( k (cid:107) , z ) ∂z − k (cid:107) (cid:98) φ ∗ ( k (cid:107) , z ) = 0 , (21)which has the analytical solution (ruling out growth at infinity) (cid:98) φ ∗ ( k (cid:107) , z ≤
0) = C e k (cid:107) z and (cid:98) φ ∗ ( k (cid:107) , z ≥ H ) = C e − k (cid:107) z , (22)where C and C are unknown constants.The form of the solution (22) implies the boundary conditions ∂ (cid:98) φ ∗ ( k (cid:107) , z → − ) ∂z − k (cid:107) (cid:98) φ ∗ ( k (cid:107) , z → − ) = 0 , ∂ (cid:98) φ ∗ ( k (cid:107) , z → H + ) ∂z + k (cid:107) (cid:98) φ ∗ ( k (cid:107) , z → H + ) = 0 . (23)Since φ ∗ and ∂φ ∗ /∂z are continuous everywhere, including at z = 0 and z = H , the boundaryconditions (23) must hold for the solution (cid:98) φ on [0 , H ] as well. These equations are the Dirichlet to Neumann map for a doubly-periodic domain [3, 4], and they give us a two point0boundary value problem (BVP) on [0 , H ] for each Fourier mode, (cid:15) (cid:32) ∂ (cid:98) φ ∗ ( k (cid:107) , z ) ∂z − k (cid:107) (cid:98) φ ∗ ( k (cid:107) , z ) (cid:33) = − (cid:98) f ( k (cid:107) , z ) , z ∈ [0 , H ] , (24) ∂ (cid:98) φ ∗ ( k (cid:107) , ∂z − k (cid:107) (cid:98) φ ∗ ( k (cid:107) ,
0) = 0 , ∂ (cid:98) φ ∗ ( k (cid:107) , H ) ∂z + k (cid:107) (cid:98) φ ∗ ( k (cid:107) , H ) = 0 . Appendix A describes the integral formulation we use to numerically solve this BVP in theChebyshev basis [16].The solution for (cid:98) φ ∗ (cid:0) k (cid:107) = , z (cid:1) is only well-defined when the integral of f ( x ) is zero.Because this is not in general the case, we will not be able to complete the formulationunless we consider both the charge densities f ( x ) and the surface charge densities σ b/t together , which we do in Section III A 3. For the moment, we consider the BVP for k (cid:107) = , (cid:15) ∂ (cid:98) φ ∗ ( , z ) ∂z = − (cid:98) f ( , z ) , (25)and define an initial solution which is only correct up to a linear mode, (cid:98) φ ∗ ( , z ) = − (cid:15) (cid:32)(cid:90) zz (cid:48) =0 (cid:90) z (cid:48) s =0 (cid:98) f (0 , s ) ds dz (cid:48) (cid:33) . (26)The linear mode will be corrected in Section III A 3.
2. Second problem: harmonic correction
We now move to the second piece of the potential in the slab, the correction φ ( c ) , thepurpose of which is to give the correct boundary conditions (7) − (10) for the total solution (cid:98) φ .Because (cid:98) φ ∗ already included the charge densities in a Poisson solve, the only charge densitiesthat remain are those on the walls. We therefore have the Laplace equation∆ φ ( c ) = 0 (27)1on −∞ < z <
0, 0 < z < H , and
H < z < ∞ , all with x and y periodic. For generality, weassume that the boundary conditions are the jump conditions φ ( c ) ( x, y, z → + ) − φ ( c ) ( x, y, z → − ) = − m ( b ) φ ( x, y ) , (28) (cid:15) ∂φ ( c ) ∂z ( x, y, z → + ) − (cid:15) b ∂φ ( c ) ∂z ( x, y, z → − ) = − m ( b ) E ( x, y ) , (29) φ ( c ) ( x, y, z → H − ) − φ ( c ) ( x, y, z → H + ) = − m ( t ) φ ( x, y ) , (30) (cid:15) ∂φ ( c ) ∂z ( x, y, z → H − ) − (cid:15) t ∂φ ( c ) ∂z ( x, y, z → H + ) = − m ( t ) E ( x, y ) . (31)Specifically, the values of m φ and m E that ensure the overall BCs (7) − (10) are satisfiedby φ = φ ∗ + φ ( c ) are given by m ( b ) φ = φ ∗ ( x, y, z → + ) − φ ∗ ( x, y, z → − ) = 0 , (32) m ( b ) E = (cid:15) ∂φ ∗ ∂z ( x, y, z → + ) − (cid:15) b ∂φ ∗ ∂z ( x, y, z → − ) + σ b ( x, y ) = ( (cid:15) − (cid:15) b ) ∂φ ∗ ∂z ( x, y,
0) + σ b ( x, y ) , (33) m ( t ) φ = φ ∗ ( x, y, z → H − ) − φ ∗ ( x, y, z → H + ) = 0 , (34) m ( t ) E = (cid:15) ∂φ ∗ ∂z ( x, y, z → H − ) − (cid:15) t ∂φ ∗ ∂z ( x, y, z → H + ) − σ t ( x, y ) = ( (cid:15) − (cid:15) t ) ∂φ ∗ ∂z ( x, y, H ) − σ t ( x, y ) . (35)These particular expressions simplify because φ ∗ is continuously differentiable across z = 0and z = H , but in Section IV we will have nonzero m φ .The solution method for the harmonic problem (27) − (31) is now straightforward. Weintroduce smooth harmonic functions φ ( c ) i , φ ( c ) t , and φ ( c ) b , whose domain is all of R , and set φ ( c ) = φ ( c ) i if z ∈ (0 , H ) φ ( c ) t if z > Hφ ( c ) b if z < . (36)2After a Fourier transform in x and y , we obtain the solution analytically as (cid:98) φ ( c ) i ( k (cid:107) , z ) = A i ( k (cid:107) ) e k (cid:107) z + B i ( k (cid:107) ) e − k (cid:107) z k (cid:107) > A i ( ) z + B i ( ) k (cid:107) = 0 , (37)with solutions of the same form for (cid:98) φ b ( k (cid:107) , z ) and (cid:98) φ t ( k (cid:107) , z ), except that boundedness of (cid:98) φ implies that A t ( k (cid:107) (cid:54) = 0) = 0 and B b ( k (cid:107) (cid:54) = 0) = 0. When k (cid:107) (cid:54) = 0, the four coefficients A i/b (cid:0) k (cid:107) (cid:1) , B i/t (cid:0) k (cid:107) (cid:1) can be determined in a straightforward way from the four boundaryconditions (28) − (31). The resulting solution in the slab interior is (cid:98) φ i ( k (cid:107) , z ) = ( (cid:15) t + 1) e − k (cid:107) z ( (cid:98) m ( b ) E ( k (cid:107) ) − (cid:15) b k (cid:107) (cid:98) m ( b ) φ ( k (cid:107) )) k (cid:107) ( (cid:15) b + 1)( (cid:15) t + 1) − ( (cid:15) b (cid:15) t + (cid:15) b + (cid:15) t − e − Hk (cid:107) (38) − ( (cid:15) b + 1) e k (cid:107) ( − H + z ) ( (cid:15) t k (cid:107) (cid:98) m ( t ) φ ( k (cid:107) ) + (cid:98) m ( t ) E ( k (cid:107) )) k (cid:107) ( (cid:15) b + 1)( (cid:15) t + 1) − ( (cid:15) b (cid:15) t + (cid:15) b + (cid:15) t − e − Hk (cid:107) + ( (cid:15) b − e ( − H − z ) k (cid:107) ( (cid:15) t k (cid:107) (cid:98) m ( t ) φ ( k (cid:107) ) + (cid:98) m ( t ) E ( k (cid:107) )) k (cid:107) ( (cid:15) b + 1)( (cid:15) t + 1) − ( (cid:15) b (cid:15) t + (cid:15) b + (cid:15) t − e − Hk (cid:107) + ( (cid:15) t − (cid:0) − e k (cid:107) ( z − H ) (cid:1) ( (cid:98) m ( b ) E ( k (cid:107) ) − (cid:15) b k (cid:107) (cid:98) m ( b ) φ ( k (cid:107) )) k (cid:107) ( (cid:15) b + 1)( (cid:15) t + 1) − ( (cid:15) b (cid:15) t + (cid:15) b + (cid:15) t − e − Hk (cid:107) The solution for k (cid:107) = 0 is only well-defined when (cid:98) φ ∗ and (cid:98) φ ( c ) are added together, as weexplain next.
3. Electroneutrality and the k (cid:107) = 0 mode For the k (cid:107) = 0 mode, (cid:98) φ ∗ ( , z ) is defined as the double integral (26). We also have linearmodes for φ ( c ) i = A i ( ) z + B i ( ) and likewise for φ ( c ) b and φ ( c ) t . The total solution (cid:98) φ ( , z ) = (cid:98) φ ∗ ( , z ) + (cid:98) φ ( c ) ( , z ) must satisfy the slab BCs (7) − (10). Substituting the representation (37)for (cid:98) φ ( c ) into the slab BCs and using (4), we obtain equations for the unknown coefficients in3terms of the double integral (cid:98) φ ∗ ( , z ), (cid:98) φ ∗ ( ,
0) + B i ( ) = B b ( ) = 0 , (39) A i ( ) H + B i ( ) + (cid:98) φ ∗ ( , H ) = A t ( ) H + B t ( ) , (40) (cid:15) (cid:32) ∂ (cid:98) φ ∗ ∂z ( ,
0) + A i ( ) (cid:33) − (cid:15) b (cid:32) ∂ (cid:98) φ ∗ ∂z ( ,
0) + A b ( ) (cid:33) + (cid:98) σ b ( ) = 0 , (41) (cid:15) (cid:32) ∂ (cid:98) φ ∗ ∂z ( , H ) + A i ( ) (cid:33) − (cid:15) t (cid:32) ∂ (cid:98) φ ∗ ∂z ( , H ) + A t ( ) (cid:33) − (cid:98) σ t ( ) = 0 . (42)There are two more equations relating to the decay of the electric field at z → ±∞ . On z < z > H , the electric field must decay to zero, which implies that the growth from φ ( c ) b and φ ( c ) t must cancel that due to (cid:98) φ ∗ on z ≤ z ≥ H , A b ( ) = − ∂ (cid:98) φ ∗ ∂z ( , z ≤ A t ( ) = − ∂ (cid:98) φ ∗ ∂z ( , z ≥ H ) . (43)Substituting this into (41) − (42), we obtain two solutions for A i ( ), which are equivalent forelectroneutral slabs, A i ( ) = − (cid:98) σ b ( ) (cid:15) − ∂ (cid:98) φ ∗ ∂z ( , z = 0) (44)= (cid:98) σ t ( ) (cid:15) − ∂ (cid:98) φ ∗ ∂z ( , z = H ) . (45)The equivalence can be seen by applying the divergence theorem to the Poisson equation(19), and then replacing the surface integrals with zero-Fourier modes to obtain ∂ (cid:98) φ ∗ ∂z ( , z = H ) − ∂ (cid:98) φ ∗ ∂z ( , z = 0) = − (cid:15) N (cid:88) k =1 q k . (46)Substituting this into the electroneutrality condition (11) yields (44)=(45). The values of B i ( ), B b ( ), and B t ( ) are determined from (39) and (40).4 B. Ewald splitting
Accurately computing the potential using the method of Section III A for a collection ofGaussian charges requires a grid spacing h ∼ g w . If g w is comparable to the distance betweencharges, then this method is sufficient to solve the problem. But as g w → without dielectric jumps orperiodicity [6, 7]. We then specialize to the doubly-periodic case with dielectric interfacesin Section IV.As discussed in the introduction, previous presentations of SE methods [5, 6, 14] viewEwald splitting as an algebraic construction for the fast calculation of electrostatic sums.The sums are split into two pieces, with the near-field sums computed by summing overneighboring points, and the far-field sums computed using non-uniform FFT methods [17].In the approach of [6], this leads to two separate Gaussian kernels: one for the algebraicsplitting, and one in the NUFFT method for spreading and interpolating to/from the FFTgrid. This is advantageous in triply periodic domains, since the NUFFT-based approachallows for the use of any NUFFT method [5, 17], including recent methods based on non-Gaussian kernels [18]. But this approach does not straightforwardly generalize to systemswith mixed periodicity. Since we instead view Ewald splitting as a smearing of charges toallow for a grid-based method, the only Gaussian kernel we use is the one in the smearingstep, and this same kernel is used both in the Ewald splitting and in the communicationbetween particles and the grid-based solver (spreading and interpolation). For triply periodicdomains, this gives an approach similar to that of [6], but it generalizes to systems withmixed periodicity straightforwardly.Beginning with the Poisson equation (6), posed in all of R , we add and subtract aconvolution (denoted by (cid:63) ) with a “screening” or “splitting” function γ . We will do the Specifically, our approach for triply periodic domains is equivalent to the SE method of [6] with theparameter η fixed to η = 1. φ = γ / (cid:63) ψ (cid:124) (cid:123)(cid:122) (cid:125) φ ( f ) + (cid:0) φ − γ / (cid:63) ψ (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) φ ( n ) , (47)where (cid:15) ∆ ψ ( x ) = − (cid:0) f (cid:63) γ / (cid:1) ( x ) . (48)The first term φ ( f ) in (47) is the “far field” and is the potential due to smeared versions ofthe original charges. We smear the Gaussian charge densities (3) using a radially symmetricGaussian kernel with standard deviation 1 / (2 ξ ) [6], γ / ( r ; ξ ) = 8 ξ (2 π ) / e − r ξ , (49)where the “splitting parameter” ξ is an arbitrary positive real number that is chosen tooptimize speed (see Section VII D).Although the far field potential defined in (47) is the convolution of γ / with ψ , we recallthat we seek the potential averaged over the Gaussian cloud g , so that our real goal is toobtain φ ( f ) = g (cid:63) φ ( f ) = (cid:0) g (cid:63) γ / (cid:1) (cid:63) ψ. (50)The kernel S := g (cid:63) γ / is a convolution of two Gaussians, one with standard deviation g w and the other with standard deviation 1 / (2 ξ ), S ( r ; g w , ξ ) = g (cid:63) γ / = 1 (cid:112) π g t exp (cid:18) − r g t (cid:19) , (51)where g t = (cid:114) ξ + g w , (52)that is, S is a Gaussian with standard deviation g t . We can now write the average far fieldpotential φ ( f ) = S (cid:63) ψ in the symmetric form φ ( f ) ( x ) = S (cid:63) ∆ − (cid:32) N (cid:88) k =1 q k S ( (cid:107) x − z k (cid:107) ) (cid:33) . (53)Here ψ ( x ) can be computed on a grid with spacing h ∼ g t using any grid-based solver, for6example an FFT-based solver for triply periodic domains.The remaining term φ ( n ) in (47) is the “near field.” Since the Laplacian commutes withthe convolution, we use (48) to write a Poisson equation for the near field (cid:15) ∆ φ ( n ) = − f + γ / (cid:63) (cid:0) f (cid:63) γ / (cid:1) = − f (cid:63) (1 − γ ) . (54)The “near field” charge density f (cid:63) (1 − γ ) for a single charge is sharply peaked, has negativetails, and integrates to zero in three dimensions. Thus the field it creates decays rapidly(exponentially) in real space, and can be computed analytically using Fourier integrals. Forthis, we recall that the charge density f for a single charge is proportional to the Gaussian g ( r ) defined in (3), with Fourier transform (cid:98) g ( k ) = exp (cid:18) − g w k (cid:19) . (55)We obtain a near field interaction kernel G ( n ) by setting f = g in (54) and solving theresulting algebraic equation in Fourier space, (cid:98) G ( n ) ( k ; g w , ξ ) = (1 − (cid:98) γ ) (cid:98) g ( k ) (cid:15)k , where (56) (cid:98) γ ( k ; ξ ) = exp (cid:0) − k / ξ (cid:1) . (57)The near field interaction kernel in real space can now be obtained by a radially-symmetricinverse Fourier transform of (56), G ( n ) ( r ; g w , ξ ) = 12 π r (cid:90) ∞ (1 − (cid:98) γ ( k )) (cid:98) g ( k ) (cid:15)k k sin ( kr ) dk = 14 π(cid:15)r (cid:32) erf (cid:18) r √ g w (cid:19) − erf (cid:32) r (cid:112) g w + ξ − (cid:33)(cid:33) . (58)In the limit as ξ → g w →
0, we recover the equation for the potential of a point chargein free space, 1 / (4 π(cid:15)r ), while in the limit g w → ξ fixed, we obtain the near fieldinteraction kernel for point charges [7]. Because the near field interaction kernel (58) decaysexponentially, it can be truncated to zero for r > r nf , where r nf ∼ g t (see Section V C 2).The near field kernel (58) will give the pointwise near field potential at any point in space7via φ ( n ) ( x ) = N (cid:88) k =1 q k G ( n ) ( (cid:107) x − z k (cid:107) ; g w , ξ ) (59)The average near field potential φ ( n ) = g (cid:63) φ ( n ) in is given by φ ( n ) ( x ) = N (cid:88) k =1 q k G ( n ) ( (cid:107) x − z k (cid:107) ; g w , ξ ) (60)where the average near field interaction kernel, G ( n ) = g (cid:63) G ( n ) . This convolution can becomputed analytically, since (cid:98) G ( n ) = (cid:98) g (cid:98) G ( n ) , to give G ( n ) ( r ; g w , ξ ) = 14 π(cid:15)r (cid:32) erf (cid:18) r g w (cid:19) − erf (cid:32) r (cid:112) g w + ξ − (cid:33)(cid:33) = G ( n ) ( r ; g w √ , ξ ) . (61)In a similar way to the potential field, we define the pointwise near field electric field by E ( n ) ( x ; g w , ξ ) = − N (cid:88) k =1 q k ∂G ( n ) ( (cid:107) r k (cid:107) ; g w , ξ ) ∂r (cid:98) r k , (62)where r k = x − z k and (cid:98) r k = r k / (cid:107) r k (cid:107) . The averaged near field electric field is given byconvolving (62) with g , or, equivalently, by differentiating the average near field Green’sfunction (61), E ( n ) ( x ; g w , ξ ) = − N (cid:88) k =1 q k ∂G ( n ) ( (cid:107) r k (cid:107) ; g w , ξ ) ∂r (cid:98) r k . (63)To avoid cancellation of digits, for small r values ( r < − g w in double precision) we usethe Taylor series for G ( n ) in (63). IV. EWALD SPLITTING FOR SLABS
In this section we develop a solver for point-like charges in a slab geometry. The firststep is to define the solution using images rather than boundary conditions at the slabwalls. For a single dielectric boundary, reflecting the charge locations across the boundaryand giving them modified strengths yields a set of image charges. The potential from theseimages plus the original charges is then what we seek. We show in Section IV A that thismethod of images can be applied to the near field and far field problems separately, with the8result that both fields satisfy the boundary conditions at the dielectric interface. Thus fora single dielectric jump, we could make a set of images and solve the problem with uniformpermittivity using Ewald splitting as described in Section III B.The situation for multiple dielectric jumps is more complex because there are infinitelymany images. The novelty in our algorithm is its ability to handle the infinitely manyimages in an efficient way. Beginning with the system of infinitely many images, we followthe Ewald splitting method of Section III B to form a near field and far field problem withthe infinitely many image charges. In the near field problem, the interaction kernel (61)decays exponentially in real space, and therefore in Section IV B we truncate it so that thecharges inside the slab only interact with their nearest images on either side of the slab.In the far field problem, the smearing of the charges and images allows for the use of agrid-based solver, and in Section IV C we evoke the components developed in Section III A tosplit the far field problem (48) into two solves. In the first solve, described in Section IV C 2,we construct an initial solution ψ ∗ for the potential inside and outside the slab. We includein this initial solve a requisite number of images so that the remaining field ψ ( c ) = ψ − ψ ∗ isharmonic and can be determined analytically in ( k (cid:107) , z ) space as described in Section IV C 3.As in Section III A, the k (cid:107) = 0 mode is only well-defined for the total potential ψ = ψ ∗ + ψ ( c ) ,and so we consider it at the end in Section IV C 4. A. Image construction
We begin by reviewing the classical image construction for point charges near a singledielectric interface at z = 0, where for z > (cid:15) and for z < (cid:15) b . Let a charge of strength q be centered at ( x, y, z > z > (cid:15) and containedthe original charge along with an image charge of strength q ∗ = − q (cid:15) b − (cid:15)(cid:15) b + (cid:15) , (64)9positioned at ( x, y, − z ). For z <
0, the potential and electric field are again as if the mediumhad uniform permittivity (cid:15) and contained a single charge positioned at ( x, y, z ), with strength q ∗∗ = q (cid:15)(cid:15) b + (cid:15) . (65)This image construction ensures that the pointwise potential φ and the normal componentof the electric displacement (cid:15) E are continuous at the wall, i.e., that the boundary conditions(7) and (9) are satisfied with σ b = 0.The same image construction can be used when the charges are spherically-symmetricclouds instead of points. To show this, let G ( r ) be the function that gives the pointwisepotential a distance r from the center of the charge. We assume that the charge is centeredat height ∆ z above the interface. Let ( x, y,
0) be any point on the xy plane, and let r be thedistance from the center of the charge to that point. Then, using the image construction,the potential at z = 0 from z < φ ( z → − ) = 2 (cid:15)(cid:15) b + (cid:15) G ( r ) , (66)while using (64) for the potential on z > φ ( z → + ) = G ( r ) − (cid:15) b − (cid:15)(cid:15) b + (cid:15) G ( r ) = φ ( z → − ) . (67)The electric displacement calculation is similar. From (65), we have (cid:15) b ∂φ∂z ( z → − ) = 2 (cid:15)(cid:15) b (cid:15) b + (cid:15) ∂G∂r ( r ) ∆ zr , (68)while using (64) we obtain (cid:15) ∂φ∂z ( z → + ) = (cid:15) ∂G∂r ( r ) (cid:18) ∆ zr − (cid:18) (cid:15) b − (cid:15)(cid:15) b + (cid:15) (cid:19) ( − ∆ z ) r (cid:19) = (cid:15) b ∂φ∂z ( z → − ) , (69)which confirms that the potential and electric displacement are continuous across the in-terface. Since the form of the kernel G ( r ) does not matter, this implies that the classicalimage construction applies to any isotropic charge clouds (in particular, Gaussian clouds),even ones that overlap the interface.0 B. Near field in the slab geometry
For the doubly periodic slab geometry, the rapid decay of the near field kernel can be usedto simplify the number of required images in the near field problem. Let us consider firstthe case of a single wall to simplify the argument. Recall that r nf is the distance at whichthe near field kernel is truncated. Then in order for the near field pointwise potential φ ( n ) to satisfy the boundary conditions at the wall, we must include all images that are centered r nf or closer to the wall.The case of the slab is similar. In order for the near field pointwise potential φ ( n ) to satisfyboundary conditions at a wall, we only need to include the images (with respect to that wall)that are r nf or closer to it. Suppose that we want to limit the near field computation tothe first set of images above/below each wall (we analyze the implications of this choice inSection V C 2). Then, as shown in Fig. 1, the key requirement is that images below thebottom wall cannot affect the top wall. Now let h be the minimum distance between thebottom wall of the slab and the center of any charge. Then the closest an image can be tothe bottom wall is h , and we have the requirement that r nf < H + h for only the first setof images below z < z > H and the bottom wall. We also ensure that r nf < min ( L x / , L y / n img gives the requirement r nf < n img H + h . However, we will show in Section V C that far-field constraints on theEwald parameter ξ give an algorithm for which r nf < H + h automatically, and so we willset n img = 1. C. Far field in the slab geometry
We are now ready to tackle the far field problem. For the dielectric slab, we need tosolve the far field Poisson equation (48) with periodicity in the x and y directions and slab1 Figure 1:
Images of the original charges in the slab (blue circles) in the near field problem.If we want the jump BCs (7) − (10) to be satisfied with σ t = σ b = 0 and only one set of imagesabove/below the slab, then the near field created by the images below z = 0 (filled red circles)cannot extend above z ≥ H , and similarly for the images above z > H (empty red circles)and the bottom wall.boundary conditions (7) − (10). To do this, we define ψ piecewise by ψ ( x, y, z ) = ψ i ( x, y, z ) 0 ≤ z ≤ Hψ b ( x, y, z ) z < ψ t ( x, y, z ) z > H , (70)where ψ i , ψ b , and ψ t are smooth fields defined for all z ∈ R . Specifically, the potentials ψ b and ψ t are due to the far-field charges and their (infinitely many) images above z = H and below z = 0, respectively. Using (cid:101) f := f ∗ γ / to denote the charge density from the(smeared) far-field charges, (cid:101) f (img) z>H to denote the density from the smeared images centeredabove the top wall, and (cid:101) f (img) z< to denote the density from the smeared images centered belowthe bottom wall, the potentials ψ b and ψ t satisfy the Poisson equations (cid:15) ∆ ψ b = − (cid:15)(cid:15) b + (cid:15) (cid:16) (cid:101) f + (cid:101) f (img) z>H (cid:17) , (71) (cid:15) ∆ ψ t = − (cid:15)(cid:15) t + (cid:15) (cid:16) (cid:101) f + (cid:101) f (img) z< (cid:17) , (72)with periodic BCs in the x and y directions and free space BCs in the z direction (i.e., nodielectric boundaries and decay of the potential as z → ±∞ ).2The interior potential ψ i satisfies the Poisson equation (cid:15) ∆ ψ i = − (cid:16) (cid:101) f + (cid:101) f (img) (cid:17) , (73)with periodicity in the x and y directions. Here (cid:101) f (img) = (cid:101) f (img) z>H + (cid:101) f (img) z< is the charge densitydue to the (infinite number of) images above and below the slab. Note that the strength ofthese images decays to zero as z → ∞ , so the Poisson equations (71) − (73) are well-posed.If we use free space boundary conditions in the z direction for (73), then we know that ψ asdefined in (70) will satisfy the jump BCs (7) − (10) with σ b = σ t = 0 (i.e., with no wall-boundcharge).Assuming that the densities σ b/t are smooth enough to be resolved by the grid used forthe far field solve, it makes sense to include the potential from the wall-bound densities,which we denote by φ ( c ) , as part of the far field potential φ ( f ) . On the slab interior, wedefine the far field potential as φ ( f ) := γ / ∗ ψ i + φ ( c ) i , (74)where φ ( c ) i is the harmonic potential defined in (36) and (37) with m φ = 0 and m E = ± σ b/t .Because φ ( c ) i is harmonic, by the generalized mean value theorem it is indistinguishablefrom its convolution with the radially-symmetric kernel γ / . This means that the far fieldpotential can be rewritten as φ ( f ) = γ / ∗ (cid:16) ψ i + φ ( c ) i (cid:17) , and φ ( f ) = S ∗ (cid:16) ψ i + φ ( c ) i (cid:17) . (75)For convenience of notation, we redefine ψ i to refer to the combined potential ψ i + ψ ( c ) i fromthe smeared charges and wall-bound densities. Using this definition, we obtain boundary3conditions for (73), ψ i ( x, y, z = 0) − ψ b ( x, y, z = 0) = 0 , (76) (cid:15) ∂ψ i ∂z ( x, y, z = 0) − (cid:15) b ∂ψ b ∂z ( x, y, z = 0) = − σ b ( x, y ) , (77) ψ i ( x, y, z = H ) − ψ t ( x, y, z = H ) = 0 , (78) (cid:15) ∂ψ i ∂z ( x, y, z = H ) − (cid:15) t ∂ψ t ∂z ( x, y, z = H ) = σ t ( x, y ) . (79)
1. Potential and charge splitting
While the Poisson equations (71) − (73) are periodic in x and y and unbounded in z , wecannot directly apply the techniques of Section III A to solve them because the image chargedensities in (71) − (73) are unbounded in the z direction. To work around this, we split eachpotential ψ i/t/b = ψ ∗ i/t/b + ψ ( c ) i/t/b (80)into the potential due to a finite number of images close to the slab (denoted with ∗ andobtained by the techniques of Section III A 1), and the harmonic potential due to the restof the (infinitely many) images (denoted with ( c ) and obtained using the method of SectionIII A 2). A conceptual picture of our far field solver is depicted in Fig. 2.To more precisely quantify the domain where we need to evaluate ψ i , we define thetruncation radius of the smeared (far-field) Gaussian charges as H E , where typically H E > g t . If the charge centers are contained on z ∈ [0 , H ], then we need to compute ψ i for z ∈ [ − H E , H + H E ] so that we can average it with the kernel S to obtain ψ i = S (cid:63) ψ i . Thusour goal is to obtain ψ i for at least z ∈ [ − H E , H + H E ] (gray region in Fig. 2).We first split the smeared charge density (cid:101) f into a piece (cid:101) f > H E coming from chargescentered at least 2 H E from the slab boundaries (i.e., charges with centers in z ∈ [2 H E , H − H E ], denoted by C > H E and shown with green circles in Fig. 2), and a piece (cid:101) f < H E comingfrom charges centered a distance 2 H E or less from the slab boundaries (denoted C < H E andshown with blue circles in Fig. 2), (cid:101) f = (cid:101) f > H E + (cid:101) f < H E . (81)4 + + Figure 2:
Far field algorithm. We split the charges and images into four groups: bluecharges = C < H E , green charges = C > H E , filled red images = C (img) < H E , dashed red circles= C (img) > H E ; the group C (img) > H E contains infinitely many images, not all of which are shown.Circles denote the support of the smeared far-field charges, solid black lines denote thedielectric interfaces, solid arrows denote a region of size 2 H E (the diameter of the far-fieldcharge support), and dotted arrows denote a region of size 3 H E (the amount we mustextend the domain to enclose the solid circles). Beginning in the left column, we firstobtain an initial guess for the solution inside the slab, ψ ∗ i , using the charges (solid blueand solid green) and images (solid red) whose support overlaps the gray region from whichwe interpolate, z ∈ [ − H E , H + H E ]. In the right column, we obtain an initial guess forthe solution outside the slab, ψ ∗ o , using the charges (solid blue) that are centered 2 H E orcloser to either wall. Both Poisson problems (for ψ ∗ i and ψ ∗ o ) are solved on a doubly periodicdomain with z ∈ [ − H E , H + 3 H E ] using the method of Section III A 1. The rest of thecharges and images (dashed circles) can be coarse-grained into a harmonic correction solvefor ψ ( c ) i and ψ ( c ) t/b .The splitting on the set of images (cid:101) f (img) = (cid:101) f (img) > H E + (cid:101) f img < H E is exactly the same, with (cid:101) f < H E representing the first set of images of C < H E , denoted by C (img) < H E and shown with solid redcircles in Fig. 2, and (cid:101) f > H E representing the rest of the images (including higher-order imagesof C < H E ), denoted by C (img) > H E and shown with dashed red circles in Fig. 2.5Using this splitting, the Poisson equation (73) for the interior solution can be written as (cid:15) ∆ ψ i = − (cid:16) (cid:101) f < H E + (cid:101) f > H E + (cid:101) f (img) < H E + (cid:101) f (img) > H E (cid:17) , (82)with the boundary conditions (76) − (79). In Fig. 2, we show (some of) the infinite number ofimages that make up the density (cid:101) f (img) > H E as dashed red circles. We see that all of these imagesare centered at least 2 H E from the slab (i.e., in the region z ∈ ( ∞ , − H E ) ∪ ( H + 2 H E , ∞ )that is outside of the dashed lines in Fig. 2), and are therefore compactly supported outsideof the gray region z ∈ [ − H E , H + H E ], on which we need to compute ψ i .
2. Intermediate potentials
We define the intermediate potential ψ ∗ i as the potential from all charges inside the slab,together with the (finite number of) images C (img) < H E whose support penetrates the region z ∈ [ − H E , H + H E ]. This is the potential due to the solid circles in the left column of Fig.2 and is given by (cid:15) ∆ ψ ∗ i = − (cid:16) (cid:101) f > H E + (cid:101) f < H E + (cid:101) f (img) < H E (cid:17) , (83)with periodicity in the x and y directions and free space BCs in the z direction. This Poissonequation can be solved using the method of Section III A 1, except that the domain in the z direction has to be modified. Because we include images centered 2 H E or closer to eitherwall, and because all smeared far-field charges have a support of H E in every direction, theright hand side (r.h.s.) of (83) is supported on z ∈ [ − H E , H + 3 H E ]. Since all imageswith support overlapping the region z ∈ [ − H E , H + H E ] are included in (83), the difference ψ ( c ) i = ψ i − ψ ∗ i must be harmonic for z ∈ [ − H E , H + H E ].Similarly, we construct solutions for ψ b and ψ t by using the splitting (80) to decomposeeach field into the potential from a finite number of charges/images and a harmonic correc-tion. In order to solve the correction problems analytically, we require that the correction ψ ( c ) b be harmonic on ( −∞ , ψ ( c ) t be harmonic on [ H, ∞ ). Thismeans that the intermediate potential ψ ∗ b must include any charges whose support extendsbelow z = 0 (those centered H E or closer to the bottom wall, like the bottom-most bluecharge in the right column of Fig. 2), and the intermediate potential ψ ∗ t must include any6charges whose support extends above z = H (those centered H E or closer to the top wall).To minimize the number of doubly-periodic Poisson solves, we combine these charges intoan intermediate potential ψ ∗ o which includes any charges C < H E centered within 2 H E of thetop and bottom boundaries (blue charges in the right column of Fig. 2), (cid:15) ∆ ψ ∗ o = (cid:101) f < H E . (84)This Poisson equation can again be solved using the method of Section III A 1. Since allof the charges included are centered within the slab and have support H E in all directions,we could use the domain z ∈ [ − H E , H + H E ]. The algorithm is simpler and more efficient,however, if we use the same domain and grid as the Poisson solve (83), z ∈ [ − H E , H +3 H E ].If we assign the charges C < H E the proper image strength (depending on whether we seek thesolution above or below the slab), we obtain the representation of the intermediate potentials ψ ∗ b = 2 (cid:15)(cid:15) b + (cid:15) ψ ∗ o , (85) ψ ∗ t = 2 (cid:15)(cid:15) t + (cid:15) ψ ∗ o . (86)
3. Correction potentials
In our initial solve (83) for ψ ∗ i , we did not include the images C (img) > H E , which are shown asdashed red circles in Fig. 2. Subtracting (83) from the original equation (82) for ψ i , we geta Poisson equation for the correction ψ ( c ) i = ψ i − ψ ∗ i which involves these images, (cid:15) ∆ ψ ( c ) i = − (cid:101) f (img) > H E . (87)Since (cid:101) f (img) > H E is compactly supported outside of z ∈ [ − H E , H + H E ] (no overlap of the dashedempty circles with the gray region in Fig. 2), we have the harmonic problem for the interiorcorrection (cid:15) ∆ ψ ( c ) i = 0 for z ∈ [ − H E , H + H E ] . (88)For the exterior corrections, the splitting of the r.h.s. (cid:101) f in (81) and definition of ψ b in7(71) imply that the correction ψ ( c ) b solves the Poisson equation (cid:15) ∆ ψ ( c ) b = − (cid:15)(cid:15) b + (cid:15) (cid:16) (cid:101) f > H E + (cid:101) f (img) z>H (cid:17) . (89)Since the r.h.s. of (89), which comprises the dashed open circles positioned above z = 0in the right column of Fig. 2, is compactly supported far away from the bottom wall (allcharges and images are centered at least 2 H E from the bottom wall), we obtain the Laplaceequation (cid:15) ∆ ψ ( c ) b = 0 on ( −∞ , . (90)The procedure for ψ ( c ) t is identical. We observe that the empty charges and images positionedbelow z = H in the right column of Fig. 2 are far from the top wall of the slab. The correctionfield that they generate is therefore harmonic above z = H , (cid:15) ∆ ψ ( c ) t = 0 on [ H, ∞ ) . (91)We now have the three Laplace equations (88), (90), and (91) (on different domains) for ψ ( c ) i , ψ ( c ) b , and ψ ( c ) t , just as in Section III A 2.The final step is to impose boundary conditions for the correction solve such that the totalfields ψ i , ψ o , and ψ t satisfy the boundary conditions (76) − (79). If we use the mismatches m ( b ) φ ( x, y ) = ψ ∗ i ( x, y, z = 0) − ψ ∗ b ( x, y, z = 0) , (92) m ( b ) E ( x, y ) = (cid:15) ∂ψ ∗ i ∂z ( x, y, z = 0) − (cid:15) b ∂ψ ∗ b ∂z ( x, y, z = 0) + σ b ( x, y ) , (93) m ( t ) φ ( x, y ) = ψ ∗ i ( x, y, z = H ) − ψ ∗ t ( x, y, z = H ) , (94) m ( t ) E ( x, y ) = (cid:15) ∂ψ ∗ i ∂z ( x, y, z = H ) − (cid:15) t ∂ψ ∗ t ∂z ( x, y, z = H ) − σ t ( x, y ) , (95)as boundary conditions for the harmonic solve described in Section III A 2, we obtain ψ ( c ) i on [ − H E , H + H E ] and, if desired, ψ ( c ) t on z > H and ψ ( c ) b on z <
0, for all modes except k (cid:107) = 0. We could combine the corrections ψ ( c ) i/t/b with the intermediate potentials ψ ∗ i/t/b toobtain ψ everywhere, although we only need ψ i on [ − H E , H + H E ] to evaluate the energyand forces on the charges. In any case, the k (cid:107) = 0 mode must be handled separately asdiscussed in the next section.8
4. Corrections for k (cid:107) = 0 In a natural extension of the method presented in Section III A 3, we obtain the solutionfor k (cid:107) = 0 by considering the total solutions ψ i/t/b = ψ ∗ i/t/b + ψ ( c ) i/t/b , rather than consideringthe intermediate and correction potentials separately. In xy Fourier space, we have thesolutions for the k (cid:107) = mode (cid:98) ψ i ( k (cid:107) = , z ) = (cid:98) ψ ∗ i ( , z ) + A i ( ) z + B i ( ) , (96) (cid:98) ψ b ( k (cid:107) = , z ) = (cid:98) ψ ∗ b ( , z ) + A b ( ) z + B b ( ) , (97) (cid:98) ψ t ( k (cid:107) = , z ) = (cid:98) ψ ∗ t ( , z ) + A t ( ) z + B t ( ) (98)By (cid:98) ψ ∗ i/t/b ( , z ), we mean the k (cid:107) = 0 components of the solutions of the Poisson equations(83) and (84), (cid:98) ψ ∗ i ( , z ) = − (cid:15) (cid:90) z − H E dz (cid:48) (cid:90) z (cid:48) − H E dz (cid:48)(cid:48) (cid:18)(cid:98)(cid:101) f < H E ( , z (cid:48)(cid:48) ) + (cid:98)(cid:101) f > H E ( , z (cid:48)(cid:48) ) + (cid:98)(cid:101) f (img) < H E ( , z (cid:48)(cid:48) ) (cid:19) , (99) (cid:98) ψ ∗ o ( , z ) = − (cid:15) (cid:90) z − H E dz (cid:48) (cid:90) z (cid:48) − H E dz (cid:48)(cid:48) (cid:98)(cid:101) f < H E ( , z (cid:48)(cid:48) ) . (100)Equation (100) defines (cid:98) ψ ∗ o ( , z ), from which (cid:98) ψ ∗ b ( , z ) and (cid:98) ψ ∗ t ( , z ) can be obtained by (85)and (86). Note that the formulas (99) and (100) imply a set of boundary conditions for thecorresponding BVPs (83) and (84) for k (cid:107) = 0; these BCs are arbitrary since (cid:98) ψ ∗ i/o ( , z ) is onlydefined up to a linear function.Once the solutions (cid:98) ψ ∗ i/b/t ( , z ) are obtained, we can solve for the coefficients A b ( ) and A t ( ) by imposing the decay boundary condition for (cid:98) ψ b in (97) and (cid:98) ψ t in (98) to obtain A b ( ) = − ∂ (cid:98) ψ ∗ b ∂z ( , z = − H E ) , A t ( ) = − ∂ (cid:98) ψ ∗ t ∂z ( , z = H + 3 H E ) . (101)The electric displacement boundary conditions (77) and (79) now give A i ( ), (cid:15)A i ( ) = (cid:15) b (cid:32) ∂ (cid:98) ψ ∗ b ∂z ( ,
0) + A b ( ) (cid:33) − (cid:15) ∂ (cid:98) ψ ∗ i ∂z ( , − (cid:98) σ b ( ) (102)= (cid:15) t (cid:32) ∂ (cid:98) ψ ∗ t ∂z ( , H ) + A t ( ) (cid:33) − (cid:15) ∂ (cid:98) ψ ∗ i ∂z ( , H ) + (cid:98) σ t ( ) . (103)9For electroneutral slabs, the right hand sides of (102) and (103) must be the same in contin-uum. Finally, as in Section III A 3, we have the continuity conditions (76) and (78) whichgive B b ( ) = (cid:98) ψ ∗ i ( ,
0) + B i ( ) − (cid:98) ψ ∗ b ( , , (104) B t ( ) = (cid:98) ψ ∗ i ( , H ) − (cid:98) ψ ∗ t ( , H ) + ( A i ( ) − A t ( )) H + B i ( ) . (105)In Section V B 2, we describe how to obtain B i ( ) in real space according to (4). The valuesof B b ( ) and B t ( ) then follow from (104) and (105). V. NUMERICAL METHOD AND ALGORITHM
In this section, we discuss the implementation of the algorithm we described in Section IV.Here we focus on getting 3 − ξ , as we explain in Section V C. A. Cutoffs and grid spacing
The splitting parameter ξ determines the truncation distances for the near field kerneland far-field charge width. Beginning with the near field, we will take the near field cutoffdistance r nf to be the minimum value such that (see (58) and (62)) (cid:12)(cid:12)(cid:12)(cid:12) ∂G ( n ) ( r nf ; g w , ξ ) /∂r∂G ( n ) ( r nf ; g w , /∂r (cid:12)(cid:12)(cid:12)(cid:12) < δ, (106) This can be verified by writing out both sides in terms of the charge densities, substituting the solutions(99) − (101), then using the identity (cid:82) H (cid:101) f (img) < H E dz = ( (cid:15) − (cid:15) b ) / ( (cid:15) b + (cid:15) ) (cid:82) − H E (cid:101) f < H E dz for far-field imagesthat overlap the bottom wall. Using a similar identity for images overlapping the top wall and theelectroneutrality condition (11) then gives the equality (102)=(103). δ = 10 − for ∼ δ = 5 × − for ∼ n σ g w . As such, we truncate the kernels (61) and (63) for r > r cut = r nf + n σ g w .We use a Fourier-Chebyshev grid to solve smooth doubly periodic problems followingthe approach of Section III A. For the far field discretization, we use FFTs in the x and y directions, so we take the (FFT-friendly) number of grid points N x/y as an input, which wecombine with g w to calculate the total far-field charge width g t using (107) and splittingparameter ξ using (52). We consider two different uniform grid spacings in x and y , with h xy = g t / . ∼ h xy = g t / . ∼ h x = L x /N x = h y = L y /N y = h xy .Although the Gaussian (51) is technically nonzero everywhere on the grid, we truncateit at a finite number of grid cells n g for efficiency in spreading the charge density to thefar field grid (because the x and y grid is uniformly spaced, this assumption allows for thepossible optimization of fast Gaussian gridding [17] in the x and y directions). This givesthe total radius of the Gaussian spreading kernel as H E = n g h xy = n σ g t , (107)where n σ is the (non-integer) number of standard deviations over which the Gaussian issupported, and the radius of support of the Gaussian charge is H E . For ∼ n g = 10, which gives a Gaussian truncated at n σ = 5 / . ≈ . n g = 12 for n σ = 6 / . ≈ . xy Fourier grid is chosen, the number of Chebyshev grid points for the z gridcan be chosen based on accuracy considerations. Specifically, we require the spacing at themiddle of the Chebyshev z grid (where it is the coarsest) to be comparable to that of theFourier xy grid. In a Chebyshev grid, the maximum spacing in the middle of the grid isat most π/ H is the total height of the domain in a doubly1periodic solve as described in Section III A, we set the number of Chebyshev points to N z = (cid:38) πH h xy (cid:39) , (108)where the integer rounding should be chosen for FFT optimality. Since the z (Chebyshev)grid is irregularly spaced, fast Gaussian gridding is not an option, and we must simply assignthe kernel (51) a value of zero at any point farther than H E in the z direction from the chargecenter. B. Algorithm
The far field solver described in Section IV C requires solving two doubly periodic Poissonequations (83) and (84), and so we begin this section by describing an algorithm to solvethem using the method developed in Section III A 1. We then summarize the completealgorithm for slabs. A GPU implementation of our algorithm is available freely at github,see https://github.com/stochasticHydroTools/DPPoissonTests/ for instructions andexamples.
1. Boundary value solver
Our doubly periodic Poisson solver is based on transforming the charge density f into itsFourier/Chebyshev representation on the grid, and consists of the following steps:1. Compute the charge density f ( x, y, z ) on the Fourier-Chebyshev grid.2. Take the 3D FFT of the charge density to obtain the Fourier and Cheyshev coefficientson the grid. We refer to this as a fast Fourier-Chebyshev transform (FFCT). See [19,c. 8] for a description of how to obtain Chebyshev coefficients using the FFT.3. Use the Chebyshev boundary value solver [16] described in Appendix A to solve theBVPs (24) for each wave number k (cid:107) using a well-conditioned integral formulation.This gives the Fourier-Chebyshev coefficients of the potential φ . Note that this stepis trivially parallelizable since each mode k (cid:107) is handled independently of others. The complex FFTs which we employ here to transform from values on the grid to Chebyshev coefficientshave 2 N z − z direction, so rounding should be such that 2 N z −
24. For the electric field, compute derivatives of φ on the grid by Fourier differentiation in x and y (i.e., by multiplying (cid:98) φ by ik x or ik y ), or differentiating the Chebyshev seriesin z . This gives the Fourier-Chebyshev coefficients of ∇ φ on the grid.5. Perform a 3D inverse fast Fourier-Chebyshev transform (IFFCT) to obtain φ and E = −∇ φ on the ( x, y, z ) grid.Note that if only energy or only forces are required then some of these steps can be skipped.
2. Ewald splitting for slabs
We now detail our Ewald splitting algorithm for computing φ and E at the centers z k of all charges in the slab geometry. The near field algorithm is straightforward. From ξ ,we determine r nf by solving (106) and then, for each charge and dielectric boundary, weconstruct one image and compute the sum (60) and/or (63) for each charge. This sumincludes only pairs of charges closer than r cut apart, and therefore the cost of the near fieldalgorithm is linear in the number of particles.The far field algorithm is more complex and worth listing in steps, with a graphic repre-sentation given in Fig. 2. Because image charges centered up to 2 H E away from the slab canbe included in some of the solves and each smeared (far-field) charge has a support of H E in all directions, we initialize a single second-kind Chebyshev grid on z ∈ [ − H E , H + 3 H E ]and a Fourier grid on [0 , L x ] × [0 , L y ]. We then perform the steps:1. Compute H E using (107). Then separate the charges into two groups: charges C < H E with z < H E or z > H − H E (solid blue circles in Fig. 2), and the rest of the charges C > H E (solid green circles in Fig. 2). Then compute the positions and strengths of thenecessary images C (img) < H E (solid red circles in Fig. 2) using (64).2. Construct the intermediate potential ψ ∗ o outside the slab by solving the doubly periodicproblem (84) as described in Section III A for z ∈ [ − H E , H + 3 H E ]. Specifically,first spread the charge q i for i ∈ C < H E onto the grid using the kernel S defined in(51). Then apply steps 2 − (cid:98) ψ ∗ o on the Fourier-Chebyshevgrid. For simplicity, use the BVP solver to also obtain the part of the k (cid:107) = 0 mode(100) by solving the corresponding BVP (84) for k (cid:107) = 0 with homogeneous Dirichlet3BCs at z = − H E and z = H + 3 H E . The intermediate potentials (cid:98) ψ ∗ b and (cid:98) ψ ∗ t inFourier/Chebyshev space can trivially be obtained by multiplying (cid:98) ψ ∗ o by the coefficientsin (85) and (86).3. Construct the intermediate potential ψ ∗ i by solving the doubly periodic problem (83)for z ∈ [ − H E , H + 3 H E ]. To do this, spread the charges q i for i ∈ C > H E and images q j for j ∈ C (img) < H E to the grid and add the result to the spreading for C < H E alreadycomputed in step 2. Then apply steps 2 − (cid:98) ψ ∗ i on theFourier-Chebyshev grid. For simplicity, use the BVP solver to obtain the part of the k (cid:107) = 0 mode (99) by solving the corresponding k (cid:107) = 0 BVP (83) with homogeneousDirichlet BCs at z = − H E and z = H + 3 H E .4. Calculate the mismatches m φ ( k (cid:107) (cid:54) = ) and m E ( k (cid:107) (cid:54) = ) given in (92) − (95) usingChebyshev differentiation in z . Use these mismatches as boundary conditions to obtainthe harmonic corrections ψ ( c ) i , ψ ( c ) b , and ψ ( c ) t analytically as outlined in Section III A 2.Evaluate the solution for each k (cid:107) and for each Chebyshev grid point z c with − H E ≤ z c ≤ H + H E , and set ψ ( c ) i ( k (cid:107) , z ) = 0 outside of this z range to avoid over/underflowerrors. Because of the ill-conditioning of the correction solve for large k (cid:107) , set ψ ( c ) i ( k (cid:107) , z )to zero for k (cid:107) > k max = π/h xy . Finally, perform N x N y independent 1D Chebyshevtransforms (FFTs) in z to obtain the Fourier-Chebyshev representation of (cid:98) ψ ( c ) i . Notethat this step can be trivially parallelized since each mode k (cid:107) is handled independently.5. To correct the k (cid:107) = 0 mode, modify the solution for ψ ∗ i ( k (cid:107) = , z ) already obtained instep 4 by adding the linear mode A i ( ) z as discussed in Section IV C 4. In the discretesetting, the equations for A i ( ), (102) and (103), give the same value for A i ( ) to onlyabout 3 relative digits, so set A i ( ) to be the mean result of the two.6. Set the far field values in Fourier-Chebyshev space to (cid:98) ψ i = (cid:98) ψ ∗ i + (cid:98) ψ ( c ) i . Then perform a3D IFFCT to obtain (cid:98) ψ i on the grid.7. Interpolate (average) ψ i at the charge centers using the kernel S given in (51) to obtainthe average far field potential φ ( f ) at the charge centers. The transforms from z space to Chebyshev space are done so that (cid:98) ψ ( c ) i can be combined with (cid:98) ψ ∗ i directlyin Fourier-Chebyshev space, and then a single 3D IFFCT performed to obtain ψ i . An alternative but lessefficient sequence is to transform (cid:98) ψ ( c ) i directly from Fourier- z space to real space via a (parallel) series of2D IFFTs in the xy plane, then add to the result from a 3D IFFCT on (cid:98) ψ ∗ i . ψ i = ψ ∗ i + ψ ( c ) i has two parts, we obtain the electric field inside the slab by differentiatingeach. We first differentiate ψ ∗ i on the grid using Fourier multiplication and Chebyshevdifferentiation. We then differentiate ψ ( c ) i analytically in z and with Fourier multiplicationfor x and y . Adding these two pieces together gives −∇ ψ i , which we then interpolateback onto the charges to obtain the average electric field E ( f ) on each charge. Note thatanother way to obtain the electric field is to interpolate the potential ψ i with the derivativeof the interpolation kernel S , but this enlarges the width of the interpolation kernel, therebyrequiring interpolation from more grid cells to obtain the same accuracy.When the energy is needed in addition to the forces, we add a constant to φ (representedby B i ( )) to set the pointwise potential φ ( ) = 0 and remove an arbitrary shift in theenergy. To do this, we sum the pointwise near field kernel at x = using the kernel (59)(with truncation at r nf and one set of images above and below the slab) to obtain φ ( n ) ( ).For the far field, we interpolate the potential ψ i with the kernel γ / at x = to obtain thepointwise far field potential φ ( f ) ( ). We then set B i ( ) in (96) so that the total potential φ ( n ) ( ) + φ ( f ) ( ) = 0. C. Constraints on the splitting parameter
Our decision to include only the first set of images in the near field and far field problemsleads to restrictions on the Ewald parameter ξ . In both of these problems, we make theassumption that images above the top wall cannot interact with points at or below thebottom wall, and vice versa. This constrains the truncation distances r nf and H E , andtherefore ξ itself.
1. Relationship between ξ and g w In order for there to be a reduction in the far field problem grid size significant enoughto justify Ewald splitting in the first place, the smearing (48) must increase the total widthof the Gaussian charge cloud by a substantial amount. Here we take substantial to meanthat the smeared Gaussian width g t is at least twice as large as the original width g w , so5that g t ≥ g w . Using (52) to solve for ξ , we get the bound g w √ ≤ ξ ⇔ ξ ≤ g w √ . (109)
2. Near field constraints
We recall the near field constraint that images above the slab cannot interact with thebottom wall, r nf < n img H + h , where h ≥ g w is the minimum distance between a chargeand one of the dielectric boundaries. This assumption gives a constraint on ξ as well. Toestimate the bound, we set r nf = α/ξ and use the upper bound on ξ from (109) to obtainthe equation δ = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂G ( n ) ( α/ξ ; 1 / ( √ ξ ) , ξ ) /∂r∂G ( n ) ( α/ξ ; 1 / ( √ ξ ) , /∂r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (110)which is independent of ξ and can be solved numerically for α to obtain a bound n img H + h > αξ > r nf . (111)The actual value of α depends on the relationship between the smeared Gaussian width g t and the original Gaussian width g w . For δ = 10 − , α ≈ . g t /g w ≥
2, while α = 3 . g t (cid:29) g w , which gives the range α ∈ [3 . , . δ = 5 × − , we get α ∈ [2 . , .
3. Far field constraints
We next consider the far field constraints on ξ . In the far field solves, we never includeimages of images explicitly, always assuming they can be included in the correction solvebecause they are a distance H + h from the slab boundaries. For this to be the case, theimages of images cannot overlap the region z ∈ [ − H E , H + H E ], meaning they must be atleast 2 H E ≥ g t away from the slab, so that in total we have the constraint H + h > H E .To estimate H E in terms of ξ , we first estimate g t . For Ewald splitting to actually reducethe computational complexity, we recall from (109) that we require g w ≤ / ( ξ √ g t , g t = (cid:114) ξ + g w ≤ ξ √ . (112)6Tolerance δ − × − Grid size h xy g t / . g t / . n g
12 10Gaussian truncation n σ /ξ < ( n img H + h ) / [3 . , .
50] ( n img H + h ) / [2 . , . /ξ < ( H + h ) / [4 . , .
95] ( H + h ) / [4 . , . Table I:
Summary of constraints for Ewald splitting for ∼ ∼ h xy , numberof grid cells n g in the support of the far-field Gaussian, and the corresponding numberof Gaussian standard deviations n σ in the truncated support. We then give bounds on 1 /ξ from the near-field constraint (111) and far-field constraint (113). Lower ends of the intervals(weaker constraints) are for the case when g w → g t = 2 g w .In the case g w →
0, we get g t = 1 / (2 ξ ). Since images of images are at least H + h awayfrom either wall, we need2 H E = 2 n σ g t ≤ n σ ξ √ < H + h ⇒ ξ < √ H + h )2 n σ , (113)or 1 /ξ < ( H + h ) /n σ in the case g w → n img = 1), the far field and near field constraint on ξ are almost identical, with the far field constraint being more restrictive. Thus the far fieldcalculation restricts the Ewald parameter so much that including more images in the nearfield is not necessary. VI. NUMERICAL TESTS
In this section, we validate our algorithm for dielectric slabs. We begin in Section VI Awith a simple verification: four charges inside of a slab unbounded in the lateral directions.In Section VI B, we verify that our answer is independent of the Ewald splitting parameter ξ by generating a reference result without Ewald splitting ( ξ → ∞ ). We conclude in SectionVI C by verifying that the force (17) is indeed the gradient of the energy (13). Throughoutthis section, we use a tolerance δ = 10 − (which gives the finer grid setting n g = 12, h xy = g t / .
4) for tests with 10 or fewer charges (Sections VI A and VI C), and a tolerance δ = 5 × − ( n g = 10, h xy = g t / .
2) for tests with more than 10 charges (Section VI B).7 z i q i (cid:16) E E ( z i ) − E F ( z i ) (cid:17) / mean (cid:16)(cid:13)(cid:13)(cid:13) E F ( z k ) (cid:13)(cid:13)(cid:13)(cid:17) (-0.17,-0.71,0.01) 1 (-9.2e-6, -1.9e-6, 2.4e-6)(0.44,-0.82,0.30) -1 (-9.4e-6, -1.9e-6, 3.0e-6)(-1.00,-0.63,1.00) 1 (-9.8e-6, -1.5e-6, 3.1e-6)(-0.40,-0.31,1.95) -1 (-9.5e-6, -1.6e-6, 2.5e-6) Table II:
Relative errors in the extrapolated values of the electric field from the Ewaldmethod as L x = L y → ∞ vs. the image method for free space, E E − E F , for four chargespositioned at z i with strength q i . A. Comparison with the image construction in free space
We first verify that our spectral Ewald method agrees with the free space solution as theperiodic xy length L x = L y = L → ∞ . The solution in free space can be derived via theimage construction with an infinite series of images (see e.g. [20]).We set up our test with four charges positioned randomly in the xy plane and with z locations 0.01, 0.30, 1.00, 1.95 (see Table II). We set (cid:15) = 1, (cid:15) t = 1 / (cid:15) b = 1 / H = 2, and g w = 10 − . We let ξ ≈ g t = 0 . H E = 0 . r nf = 1 . L = 28 (236 × ×
84 grid)and L = 32 (270 × ×
84 grid), extrapolate the result to L = ∞ (based on the fact thatthe finite-size correction decays like 1 /L ), and denote the resulting electric field by E E . Weshow the normalized difference E E − E F , where E F is the free space result, in Table II, whichshows that we obtain about 5 digits of accuracy for all four charges. B. Independence from splitting parameter
We next verify that our results are independent of the splitting parameter ξ . To do this,we compare the results for the spectral Ewald method with various ξ to reference resultsobtained without Ewald splitting ( ξ → ∞ ). We generate the reference results to 7 digits ofaccuracy using the method of Section III A with a fine grid.We consider 100 random charges of unit strength and alternating sign positioned on[ − , × [ − , × [4 . g w , H − . g w ] and set g w = 0 .
025 so that it is not too expensive togenerate the reference solution to high accuracy. For the jump conditions, we set (cid:15) = 1, (cid:15) t = 1 /
50, and (cid:15) b = 1 /
20. We take H = 0 .
75, so that, for some choices of the Ewald8 ξ Grid size H E r nf Error std. ∞ (ref.) 128 × × × ×
59 0.5 0.71 2.7e-59.2 40 × ×
71 0.25 0.35 5.1e-512.2 50 × ×
77 0.20 0.27 6.0e-526.0 76 × ×
92 0.13 0.16 7.9e-5
Table III:
Splitting parameters for Section VI B. Grid is on [ − L, L ] × [ − L, L ] × [ − H E , H +3 H E ] with H = 0 .
75 and L = 2. We give the standard deviation of the errors shown in Fig.3.parameter, charges in the center of the slab have (smeared far-field) support that intersectsboth z = 0 and z = H .We choose a set of Ewald parameters in the range ξ = 4 to ξ = 26. The lower bound on ξ is chosen so that each point interacts only with its nearest image to 4 digits in the near fieldcalculation (this includes periodic images as well as images in the z direction). The upperbound on ξ is the point where the grid is almost as large as that used for the calculationwithout Ewald splitting (i.e., ξ is so large that the Ewald splitting is useless). Table IIIsummarizes the values of the Ewald parameter we use and the resulting values of H E and r nf .For each Ewald parameter, we solve for the averaged potential and electric field andcompute the error in each component of E ( z i ) (relative to the ξ → ∞ reference result),normalized by the mean magnitude of E . We do this 10 times, so that there are a totalof 3000 observations E . Figure 3 shows histograms of the relative errors for each of theEwald parameters. For all ξ , we obtain distributions centered approximately around 0 with4 digits of accuracy in E . Table III shows that the standard deviation of the relative errorsin each case is ≈ × − . This demonstrates that the accuracy is better than the targeted δ = 5 × − , especially for smaller Ewald parameters when the near field dominates. C. Energy and force
Finally, we verify that the forces (17) are the gradient of the energy (13), even in thepresence of charged walls. To do this, we set (cid:15) = 1, (cid:15) t = 1 / (cid:15) b = 1 /
20, and H = 1 and L = 2. We place 10 charges with unit strengths of alternating sign randomly and uniformlyon [ − L, L ] × [ − L, L ] × [4 . g w , H − . g w ], where g w will be varied. We also place a Gaussian9 -3 -2 -1 0 1 2 310 -4 -3 -2 -1 0 1 2 310 -4 -4 -4 Figure 3:
Histogram of relative errors in the components of the forces on 100 randomlypositioned charges of unit strength. We compare four different values of the splitting param-eter ξ in the spectral Ewald method to a reference solution computed to 7 digits of accuracyusing the method of Section III A. The errors are normalized by the mean magnitude of theforce over all charges, and the experiment is repeated 10 times.charge density on the top and bottom walls, σ b/t ( x, y ) = (+ / − ) 12 √ π s exp (cid:18) − x + y s (cid:19) . (114)We use s = 0 . ξ = 6 . g t ≈ . H E = 0 . r nf = 0 .
48) for a grid of size 38 × ×
87 points. Sincethe grid spacing is chosen to resolve g t = 0 .
07, and the wall charge density is a Gaussianwith standard deviation s > g t , we expect that our grid will easily be able to resolve thewall density, and therefore the integrals in the second term of the energy (13). We haveconfirmed that the potential computed by the correction solve matches the analytical resultfor a potential of a Gaussian surface charge (not shown).0 g w − − − − | W − W | /W Table IV:
Relative differences in computed rates of work | W − W | /W using (115) and(116) for δ = 10 − with M = 10 charges.We fix δ = 10 − and compute a rate of work by finite differences of the energy (13) W = − U ( X + δ δ X ) − U ( X − δ δ X ) δ , (115)where δ X is a random displacement vector of unit length for each charge. We compare thisto the rate of work done by the force on the charges (17), W = N (cid:88) k =1 F k · δ X k = N (cid:88) k =1 q k E ( z k ) · δ X k . (116)To avoid cancellation of digits as g w →
0, we subtract the self potential, ¯ φ ( z i ) → ¯ φ ( z i ) − q i G ( n ) ( r = 0; g w , ξ = 0), so that the potentials do not become singular as g w →
0. As shownin Table IV, W = W to at least 3 digits even for g w = 10 − . VII. APPLICATIONS
In this section we apply our method to study the influence of polarization effects (imagecharges) on the structure of electric double layers in electrolyte solutions confined by oneor two dielectric boundaries. This is a well-studied problem and it is not our intention toprovide a detailed account of the topic. Rather, our main goal is to test our method bycomparing to published results, and to demonstrate the importance of accounting for jumpsin the dielectric permittivity. We also study the performance of the algorithm on GraphicalProcessing Units (GPUs) for realistic electrolyte parameters.We use Brownian Dynamics (BD) without hydrodynamic interactions to equilibrate anelectrolyte solution in a slit channel; this requires evaluating the forces on the charges ateach time step, but does not require the electrostatic energy. In a number of previous stud-ies of electrolytes, Markov Chain Monte Carlo (MCMC) is used for equilibration. However,MCMC with local moves is expensive for electrolytes since computing the long-ranged elec-trostatics is a global operation and has to be done after each trial move. We have not been1successful with controlling the rejection rate in MCMC methods based on global moves, in-cluding ones combining BD with MCMC like the Metropolis-adjusted Langevin Algorithm(MALA) [21, 22]. As a compromise, we utilize here a BD method specifically designed forcomputing equilibrium distributions to second-order accuracy [23]. It is important to em-phasize that because we use BD, we operate in the canonical
N V T ensemble and not in thegrand canonical ensemble. Furthermore, by construction our method requires maintainingstrict electroneutrality in every configuration. Although some conventions for handling non-neutral slabs have been proposed [24], and can straightforwardly be incorporated into ourhandling of the k x = k y = 0 mode, a unique well-defined meaning to the doubly-periodicPoisson equation cannot be given without strict electroneutrality.We begin by summarizing our Brownian Dynamics method in Section VII A. In SectionVII B we consider a monovalent binary electrolyte and study the depletion of charges neara water-air interface due to repulsion by the image charges. In Section VII C we study thesame effect but in a slit-channel with charged dielectric walls and with only the counterionpresent. Finally, in Section VII D we discuss the performance of our algorithm and theoptimal choice for the Ewald splitting parameter. A. Brownian Dynamics simulations
We consider a solution of monovalent ions with effective radius a and charge ± e solvatedin water. Unless stated otherwise, simulations are carried out at room temperature ( T =298 K). We include a steric repulsion between ions in the form of a modified, truncated, andregularized Lennard-Jones (LJ) potential [8] with repulsive force F steric ( r ) = − ∂U steric /∂r , F steric ( r ) = F LJ ( r = r m ) r ≤ r m F LJ ( r ) r m < r ≤ /p (2 a )0 r > /p (2 a ) (117)where r is the ion-ion distance and F LJ = − ∂U LJ /∂r where U LJ ( r ) = 4 U (cid:32)(cid:18) ar (cid:19) p − (cid:18) ar (cid:19) p (cid:33) + U . (118)2Here r m denotes a cutoff distance below which we molify the divergence of the classicalLJ potential to avoid very large forces which lead to numerical instabilities [8]. In somesimulations, we also decrease the exponent p from the traditional p = 6, in order to makethe potential softer and thus increase the stable time step size. The value of the stericrepulsion U is chosen to be ∼ k B T to avoid overlap but also avoid very large forces uponoverlap. We add a steric repulsion with walls (dielectric boundaries) by considering a virtualimage charge and computing the steric repulsion with that image charge.We employ the lower precision (3 digits of accuracy) set of parameters for our algorithm(see last column in Table I); we have confirmed that the forces on the charges did notchange to at least 3 digits upon changing the Ewald splitting paramater ξ . The algorithmis implemented on GPUs using the Universally Adaptable Multiscale Molecular Dynamics(UAMMD) [25] framework, and we use single precision in all calculations, except that thecorrection potential (38) is computed in double precision to avoid overflow and underflow.We have confirmed that switching all calculations to double precision does not lead to adegradation of the overall accuracy of the method. Our GPU code is freely available onthe UAMMD github https://github.com/RaulPPelaez/UAMMD , and the tests given in thissection can be reproduced using the codes and input files available at https://github.com/stochasticHydroTools/DPPoissonTests/ .We solve the stochastic equations of Brownian Dynamics without hydrodynamic interac-tions, d z = µ ( F ( z ) + F steric ( z )) dt + (cid:112) k B T µ d B (119)where the electrostatic forces F are defined in (17), µ is a scalar mobility, and B ( t ) isa collection of independent Brownian motions (Wiener processes). These equations areintegrated in time using a non-Markov variation of the classical Euler-Maruyama methodproposed in [26] and studied theoretically in [23]. The scheme updates the ion positionsfrom time step n to n + 1 using z n +1 = z n + µ ( F steric ( z n ) + F ( z n )) ∆ t + (cid:114) k B T µ ∆ t (cid:0) W n + W n +1 (cid:1) , (120)where W n is a collection of independent Gaussian random variables of mean zero and vari-ance one. The scheme (120) has first order weak accuracy at finite times but transitions3to second-order weak accuracy exponentially fast in time [23]. In particular, it approx-imates expectation values with respect to the Gibbs-Boltzmann equilibrium distribution ∼ exp ( − ( U + U steric ) / ( k B T )) to second-order accuracy in the time step size ∆ t . We havefound this scheme to be the most accurate among several we tested, and also found it to bejust as robust for larger time step sizes as the less-accurate Euler-Maruyama scheme. Wereport the values of ∆ t and time in units of the diffusive time τ = a / ( k B T µ ), so thatthe precise value of µ is irrelevant. Charges are checked every step to be inside the alloweddomain (between z = 0 to z = H ). If a charge leaves the domain after a given updatethe step is discarded and retried until it results in a valid configuration. After a maximumallowed number of retries an unrecoverable configuration is declared to avoid stagnation.We choose ∆ t to be small enough such that steps are rejected very rarely and unrecoverableconfigurations do not occur over long periods of time. We find that we can use a larger timestep size if we restrict the maximum allowed displacement per time step for a particle to amultiple of the ion radius; we return to this point in the Conclusions.It is straightforward to adapt our method and implementation from a doubly-periodic toa triply-periodic domain. In Fig. 4 we show results for the pair correlation function g ( r ) ofcounter ions in a dilute binary monovalent electrolyte of molarity 0.05M in a triply periodicdomain. The parameters for the simulations are given in Table V and in the caption ofthe figure. In particular, we study two different systems in the remainder of this section.The first system (second column of table) is meant to correspond to the hard-sphere systemstudied theoretically and using Monte Carlo simulations in [27]. To mimic hard spheres,we set the exponent to p = 6 in the LJ potential (118) and choose U to be sufficientlylarge to make g ( r ) almost zero for r < a . The second system (third column of table)we study uses a softer repulsion with exponent p = 2 to avoid stiffness. In this secondsystem, however, we wish to mimic point charges closely so we set g w (cid:28) a , which leads tostiffness of the electrostatic potential. When neither strict steric exclusion or point-chargeelectrostatics is required, one can use a small p = 2, large r m ∼ a , and large g w ∼ a , tomake the combined interaction potential non-stiff and allow for a larger ∆ t ; we return tothis point in the Conclusions.For comparison to the numerical results, in Fig. 4 we also show theoretical predictionsfor the pair correlation function based on Debye-Huckel-Onsager theory, modified to account4Figure 5 6 a . L xy . a aH a ag w . a . aU . k B T . k B Tr m . a ap Table V:
Parameters used in the BD simulations corresponding to Figs. 5 and 6. The timestep size is ∆ t = 2 · − τ .for the fact that our charges are not point charges, g DHO ( r ) = exp (cid:18) − U steric ( r ) + U ( r ) exp( − r/λ ) k B T (cid:19) (121)where U ( r ) is the electrostatic potential between two counter ions with Gaussian chargedensity (see first term in (61)), U ( r ) = − e π(cid:15)r erf (cid:18) r g w (cid:19) . (122)Here the Debye length is λ = (cid:18) (cid:15)k B Te N a (2 M ) (cid:19) , (123)where N a is the Avogadro number and M is the salt molar concentration (molarity) inmol /m . The numerical results in Fig. 4 are in good agreement with the theory especiallyin the Debye (exponential) tail (see inset). For the system with p = 6 we also show thetheoretical g ( r ) if the steric repulsion were a hard-sphere potential (strictly no overlaps). B. Uncharged wall
In this subsection we reproduce results from the article [27] using the doubly-periodicPoisson solver in UAMMD. We simulate an electroneutral monovalent solution inside a slabwith uncharged walls and different permittivities above and below it; the parameters arelisted in the second column in table V. We compare our results and those obtained in [27] inFig. 5. Since the walls are uncharged there is a symmetry between the anions and cations,5 -2 -1 g (r)- r / a g ( r ) UAMMD p=2p=2 DHOUAMMD p=6p=6 DHOHS DHO
Figure 4:
Pair correlation function (symbols) of counter ions in a dilute binary monovalentelectrolyte of molarity 0.05M in a triply periodic cubic domain ( N = 8192 ions), for twodifferent systems ( p = 6 or p = 2 in the LJ potential (118)). The parameters for the potential( p , r m , U , and g w ) and time step size ∆ t are given in Table V. Solid lines are the DHOtheory (121), and the dashed line shows the theory for a hard-sphere steric repulsion withthe rest of the parameters as for the system with p = 6. The inset shows the exponentialtail ∼ exp( − r/λ ) /r in g ( r ) − n ( z ) averaged among the two types of ions.There are a number of recent theories that account for the polarization effects next toa dielectric boundary, but all of them involve nontrivial computations — for recent exam-ples see [28] without and [29] with steric repulsion. The authors of [27] developed severaltheoretical approaches to compute the equilibrium density of ions next to a planar interface(wall) with hard-sphere steric repulsion in addition to electrostatics (i.e., the so-called prim-itive model of electrolytes), with the “BBGKY+EN” theory being in best agreement for themolarities we study here. We have extracted the theoretical curves from the figures in [27]and show those for comparison in Fig. 5.The MCMC simulations and theory in [27] used the grand canonical ensemble and con-6 z / a n ( z ) / n b UAMMD 1.0MTheory 1.0MUAMMD 0.5MTheory 0.5M
Figure 5:
Number density of monovalent ions n ( z ) at a distance z from a dielectric bound-ary, normalized by the “bulk” density far away from the wall. The results from our codeUAMMD are compared with the results in figures 3 and 4 in [27]. The bulk ionic density is n b = 0 .
495 mol/L = 2 . · − a − (black line) or n b = 1 .
00 mol/L = 5 . · − a − (orangeline); other parameters are listed in the second column in Table V. The total number ofions is N = 9398 (red dots) and N = 19142 (blue dots). Simulations were initialized byrandomly placing ions in the region a < z < H − a and equilibrated for 20 diffusive times τ , and then data was collected every diffusion time for the subsequent 200 diffusive times.Error bars (comparable to or smaller than the symbol size) come from averaging over 32independent simulations.sidered a system that has a bottom wall but is unbounded in the z direction, i.e., a systemthat is in contact with an infinite reservoir as z/a → ∞ with bulk molarity M bulk . Sincewe use the canonical ensemble, we put a top wall at z = H = 50 a but set (cid:15) t = (cid:15) in orderto mimic an “open” reservoir boundary. We numerically integrated the theoretical curves in[27] to obtain a total number of ions that would give the same bulk molarity far away fromthe wall, assuming that the density was constant for z (cid:29) a (which is in agreement withsimulation results).7The results in Fig. 5 show excellent agreement between our numerical results and thetheoretical predictions from [27], both for bulk molarity M bulk ≈ .
5M and M bulk ≈ M bulk = 1M predicted by the theory, and confirmed by the simulations.In this example, the depletion layer extends only over a distance of a couple of ion diametersaway from the wall. C. Charged wall test
In this subsection we study a slit channel with charged walls with only the counterionspresent in the interior of the channel. Specifically, we place N ions of charge + e in thechannel and keep the system electroneutral by setting the surface charge densities at thetop and bottom walls to σ b = σ t = σ = − ( N e ) / (cid:0) L xy (cid:1) ; the parameters are listed inthe third column in Table V. For this special case it is possible to analytically solve thePoisson-Nernst-Planck (PNP) equations for the equilibrium number density of ions acrossthe channel (see Section 4.1 in [30]). If we denote the sterically-accessible thickness of thechannel as d = H − a and solve the PNP equations with charged walls inside a channel ofwidth d , we obtain the ion density profile for a < z < H − a (zero outside of this interval), n ( z ) = n m cos ( K ( z − H/ /d ) , (124)where the dimensionless constant K < π depends on the charge density, σ = − e (cid:90) H − az = a n ( z ) dz = − K tan ( K/ (cid:15)k B Tde , (125)and n m is the number density of ions in the middle of the channel, n m = 2 K (cid:15)k B Td e . (126)In Fig. 6 we compare numerical results for n ( z ) with the prediction of the PNP equations(124) at two different densities. We see an accumulation of ions near the boundaries becauseof the surface charge, but the maximum density depends on the dielectric contrast. We8consider three different materials on the outside of the channel: glass (permittivity ≈ (cid:15) t = (cid:15) b = (cid:15) out = (cid:15) . This is not unexpected since the PNP equation does not take into accountpolarization effects. Namely, for an electroneutral slit (or cylindrical) channel, Gauss’s lawdictates that the electric field just outside of the channel must be zero; the electric fieldscreated by the ions on the outside of the channel are evanescent and vanish in the meanfield approximation. Only if charge electroneutrality is broken does the dielectric constanton the outside of the channel begin to matter in the PNP equations [31]. While polarizationeffects can be included in modified PNP equations that take into account charge correlations[28, 29], the resulting nonlocal equations are no longer simple to solve even numerically. Itis important to note that the mismatch between the PNP theory and the numerical resultsseen in Fig. 6 extends over a relatively large distance from the wall, on the order of tensof particle radii, emphasizing the need to account for dielectric jumps in theoretical andnumerical descriptions of confined electrolytes.The results in Fig. 6 show that there is only a small difference between the results obtainedfor glass and a fictitious material with (cid:15) out = 0, for which there is no dielectric displacementoutside of the channel. This is good news for approaches such as immersed boundary meth-ods [8] that only resolve the electric fields inside the channel, and impose (inhomogeneous)Neumann conditions on the electric field, rather than the correct jump boundary conditionsfor the dielectric displacement. D. Computational performance
In this section we study the computational performance of our GPU implementation ofthe doubly-periodic (DP) method, as a function of the Ewald splitting parameter ξ . Werandomly and uniformly place 2 · charges inside a slab of extent H = 50 and length L xy = 185, and set g w = a/
4. We use single precision (except for the correction solveto avoid overflow) and time each part of the algorithm using profiling tools; this requiresthat all kernel launches are synchronous, not allowing for overlapping CPU/GPU work andconcurrent kernels (which improve performance only slightly).9 (z - H/2) / a n / n m PNP (K=2)UAMMD ε out =0UAMMD ε out =5 ε UAMMD ε out = ε (z - H/2) / a n / n m PNP (K=2.65)UAMMD ε out =0UAMMD ε out =5 ε UAMMD ε out = ε Figure 6:
Number density of cations n ( z ) inside a slab channel with oppositely chargedwalls, at two different surface change densities (lower in the left panel). The density isnormalized by the density in the middle of the channel predicted by the PNP equations (see(126)), n m = 3 . · − a − (left, N = 6140 ions) and n m = 1 . · − a − (right, N = 1739ions); other parameters are listed in the third column in Table V. Numerical results obtainedusing the UAMMD doubly-periodic code are shown for three different permittivities outsideof the channel (glass, water, and an unphysical zero permittivity medium). To accelerateequilibration, initial configurations were generated by randomly and uniformly placing ionsand using rejection to make the density vary with z according to (126). We equilibrated for20 diffusive times τ , and then data was collected every diffusion time for the subsequent 200diffusive times. Results are averaged over 64 independent runs; the error bars are smallerthan the symbol size.The total wall clock time to compute the force on and potential at each charge, and thebreakup among different components of the algorithm, are shown in the left panel of Fig.7. We isolate the time due to the following components of the algorithm: Computing thecorrection field to account for the missing images, spreading the charges to the grid, com-puting the near-field interactions, interpolating the field back onto the charges, computingthe Fourier-Chebyshev Transform (FCT) using 3D FFTs, and solving the boundary valueproblems (BVPs) in the doubly-periodic Poisson solver.For comparison, in the right panel of Fig. 7 we also give timing results for a triply periodic(TP) domain (also available in UAMMD [25]), for which the Poisson equation can be solvedentirely in Fourier space using 3D FFTs. We use the same particle configurations as forDP domains in the TP domain, but we set the length of the domain in the z direction to L z = 2 H , which makes the sizes of the grids used in the 3D FFTs similar. Unlike the DPsolver, the TP solver does not require any images, a correction solve, or BVP solver, and,importantly, uses a uniform grid for which spreading and interpolation are faster in the z T i m e ( m s ) NxyξHBVPFCTInterpolationNearSpreadCorrection 0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 32 42 54 64 70 80 88 96 1087.5 12.5 17.5 22.5 27.5 T i m e ( m s ) Nxy2ξHConvolveFFTInterpolateNearSpread
Figure 7:
Execution time for different components of the algorithm for 2 · charges asa function of the Ewald splitting parameter ξ (top axis), which controls the grid size in the xy directions (bottom axis). The left panel is for a doubly-periodic (DP) domain, while theright one is for a triply-periodic (TP) domain. Timings are collected using the UAMMD[25] code compiled and run on a RTX2080Ti GPU using NVIDIA’s CUDA 11.0.direction. At the optimum split (grid size N xy = 88 for DP and 70 for TP), the time it takes tocomplete the calculation is 0 .
84 ms for TP and 4 . VIII. CONCLUSIONS
We have developed a spectrally-accurate fast method for computing electrostatic energyand forces for a collection of charges in doubly-periodic slabs with jumps in the dielectricpermittivity at the slab boundaries. To do this, we used a modification of the Spectral Ewald In particular, fast Gaussian gridding [17] can be used in all directions, although we do not find this tolead to substantial improvement on the GPU. On the GPU, the main savings comes from the fact thatthe number of grid points to spread to or interpolate from in the z direction is constant (10 grid cells ineach direction) in the TP case. smooth mismatch in theboundary conditions (BCs) at the dielectric interfaces using the analytic solution of theLaplace equation with inhomogeneous boundary conditions. Combining the far-field andcorrection steps yields a method that requires the same components as the Spectral Ewaldmethod for triply-periodic domains: spreading and interpolation from the charges to the gridusing Gaussian kernels (but note that the Chebyshev grid is not uniform in the z direction),and a forward and inverse three-dimensional FFT. All components can be parallelized onGPUs, and our public-domain implementation takes only ∼ xy direction and thus compute electrostatic interactions by summing over only thenearest image when computing electrostatic energy. Our results in Fig. 5 appear to be ingood agreement with the results given in [27], suggesting this approximation was justifiedfor the system studied. However, we also study here a case where there are no counterions,and there is no screening, so that it is crucial to properly implement periodicity in the xy directions.Another important approximation that greatly simplifies the problem is to assume that2the dielectric constant is zero outside of the slit channel. While this is unphysical (thesmallest possible permittivity is that of vacuum), the dielectric constant of most materialssuch as glass or lipid membranes is much lower than that of water, which is the typicalsolvent. In the unphysical limit of zero dielectric constant in the exterior of the channel,instead of jump conditions on the electric displacement we get a Neumann condition on theelectrostatic potential on the interior of the channel. This means that it is no longer necessaryto compute the electrostatic potential outside of the channel or worry about infinitely manyimages, similar to the case of metallic boundaries. Our results in Fig. 6 suggest that fora water-glass interface the approximation of zero permittivity outside the channel is quiteaccurate.It is straightforward to incorporate metallic slab boundaries (electrodes) in our methodsinstead of dielectric jumps; only minor modifications to the correction solve are required toset the potential at the top and/or bottom electrodes to a specified value. Since the case oftwo dielectric jumps is the hardest and requires the full power of our method, we focused onthis case in our tests.Our experience with BD for electrolyte solutions suggests that an important problem toovercome in future work is the small time step size required to stably integrate a system withstiff steric repulsion. To mimic hard-sphere repulsion in the presence of walls, we required∆ t ≈ · − τ , where τ is the typical time it takes an ion to diffuse a distance of its radius.Even with the algorithmic improvements we developed here and the computational power ofGPUs, this ∆ t is too small to reach time scales relevant to dynamics of electrolytes, includingrelaxation and charging dynamics of double layers, and electrohydrodynamic phenomena.As in recent work based on the Discrete Ion Stochastic Continuum Overdamped Solvent(DISCOS) method [8], here we introduced several mollifications that helped increase ∆ t .The first one was to mollify the charges by making them Gaussian clouds of finite widthinstead of point charges, which helps avoid the r − divergence of electrostatic forces forpairs of ions. We also mollified the traditional Lenard-Jones steric repulsion by capping itsdivergence at overlap and softening the potential by lowering the exponent to p <
6. Thesetwo changes are related to each other and have to be done in tandem because the stericrepulsion (mimicking Pauli exclusion) has to be strong enough to prevent overlap of counterions.The primary source of instability in temporal integration appears to be large steric or3electrostatic forces that occur upon ion overlap. These large forces occasionally lead todisplacements that are several ionic radii large, which is particularly problematic for ionsnear walls, as the displacements can lead to ions leaving the domain. We have had somesuccess increasing the time step size by limiting the largest possible displacement of an ionduring a time step to a fraction of the ionic radius. Even with this ad hoc change the largeststable and accurate time step size we achieved was ∆ t ≈ − τ , which is still quite smalleven though it is much larger than ∆ t in MD. We hope that future mathematical study oftemporal integrators for overdamped BD with stiff steric repulsion will lead to improvementsand allow us to reach a desirable ∆ t ≈ − τ . Another avenue worthy of exploration isreusing information between time steps, especially for the far-field Poisson solve. Since thefar-field potential and fields are smooth, it is perhaps possible to not repeat all steps of thesolve at each time step. Our preliminary investigations found some promise in this directionfor dilute electrolytes.Another important direction is to generalize our method to Stokes flow so that hydro-dynamic interactions can be accounted for in Brownian Dynamics. Our doubly-periodicFourier-Chebyshev solver can straightforwardly be generalized to the Stokes instead of thePoisson equation, as we will present in future publications. The main challenge is Ewaldsplitting in the presence of no-slip boundaries (bottom wall only or top and bottom walls);for triply periodic systems one can use the Positively Split Ewald method [32] or relatedSE methods [33] that rely heavily on Fourier transforms in all directions. Some progress onreal-space based Ewald splitting with boundaries has been made in [34]; however, becausean image construction was not used to handle the boundaries, the near field does not satisfythe BCs on the wall (as it did in our method for the Poisson equation). While the mis-match in BCs can in principle be fixed with a correction solve [34], the grid required for anaccurate correction would be much finer than the grid used for the far-field solver, negatingthe advantages of Ewald splitting. Furthermore, the method in [34] cannot handle a singlebottom wall as does our approach based on the Dirichlet to Neumann map.Ewald methods for Stokes flow based on image constructions have been developed for asingle bottom wall using Fourier transforms in all directions [35] or FMMs [36, 37]. However,the image construction for a no-slip wall for Stokes flow involves several types of imagesingularities [38], and this leads to substantial complexity and inefficiency compared to theapproach we developed here for the Poisson equation. It should be mentioned that recent4investigations using the DISCOS method demonstrate that hydrodynamic interactions makeimportant contributions to transport in electrolytes [8]. That said, these recent studies alsodemonstrate that hydrodynamic interactions can be coarse grained at scales smaller than thetypical ion-ion distance and replaced by standard non-hydrodynamic or “dry” diffusion (as weused in this work), which suggests that Ewald splitting may not be necessary for electrolytesolutions. Nevertheless, it remains a challenge for the future to adapt the method developedhere to Stokes flow. ACKNOWLEDGMENTS
We thank Zecheng Gan for helpful discussions regarding electrostatic energy in the pres-ence of surface charges. Ondrej Maxian is supported by the National Science Foundation(NSF) via GRFP/DGE-1342536. This work was also supported by the NSF under awardDMS-2011544 and through a Research and Training Group in Modeling and Simulationunder award RTG/DMS-1646339. Ra´ul P. Pel´aez acknowledges funding from Spanish gov-ernment MINECO project FIS2017-86007-C3-1, and thanks Prof. Rafael Delgado-Buscalionifor his support and additional funding.
Appendix A: Boundary value solver
Our boundary value problem (BVP) solver for (24) is based on the specgtral integrationmethod of [16]. Without loss of generality, it is most convenient when reviewing this for-mulation to assume the domain is z ∈ [ − , y (cid:48)(cid:48) ( z ) − k y ( z ) = f ( z ) , (A1) y (cid:48) (1) + ky (1) = α, y (cid:48) ( − − ky ( −
1) = β. (A2)An analytical solution can be derived for this BVP, but it requires numerically computingan integral with integrand related to e kz . Since this calculation must be done using very finegrids for kz (cid:29)
1, we prefer to use a well-conditioned more general BVP solver.5Now, let us expand all functions in truncated Chebyshev series y ( z ) = N − (cid:88) n =0 (cid:98) y n , y (cid:48) ( z ) = N − (cid:88) n =0 (cid:98) y (cid:48) n T n ( z ) , y (cid:48)(cid:48) ( z ) = N − (cid:88) n =0 (cid:98) y (cid:48)(cid:48) n T n ( z ) , f ( z ) = N − (cid:88) n =0 (cid:98) f n T n ( z ) . (A3)To obtain the coefficients (cid:98) f n , we can evaluate f on a Chebyshev grid with N points, take theperiodic extension of f ( z ), and do a complex FFT [19]. It can be shown, using the indefiniteintegrals of Chebyshev polynomials, (cid:90) T ( x ) dx = T ( x ) + C, (A4) (cid:90) T n ( x ) dx = 12 (cid:18) T n +1 ( x ) n + 1 − T n − ( x ) n − (cid:19) + C, n ≥ , (A5)where C is constant, that (cid:98) y (cid:48) = 12 (2 (cid:98) y (cid:48)(cid:48) − (cid:98) y (cid:48)(cid:48) ) (cid:98) y (cid:48) n = 12 n (cid:0)(cid:98) y (cid:48)(cid:48) n − − (cid:98) y (cid:48)(cid:48) n +1 (cid:1) , n ≥ . (A6)Computing the coefficients of the function y ( z ), we have the relationships (cid:98) y = 12 (2 (cid:98) y (cid:48) − (cid:98) y (cid:48) ) = (cid:98) y (cid:48) −
18 ( (cid:98) y (cid:48)(cid:48) − (cid:98) y (cid:48)(cid:48) ) , (cid:98) y = 14 ( (cid:98) y (cid:48) − (cid:98) y (cid:48) ) = 14 (cid:18)
12 (2 (cid:98) y (cid:48)(cid:48) − (cid:98) y (cid:48)(cid:48) ) −
16 ( (cid:98) y (cid:48)(cid:48) − (cid:98) y (cid:48)(cid:48) ) (cid:19) , (cid:98) y n = 12 n (cid:0)(cid:98) y (cid:48) n − − (cid:98) y (cid:48) n +1 (cid:1) = 12 n (cid:18) n − (cid:0)(cid:98) y (cid:48)(cid:48) n − − (cid:98) y (cid:48)(cid:48) n (cid:1) − n + 2 (cid:0)(cid:98) y (cid:48)(cid:48) n − (cid:98) y (cid:48)(cid:48) n +2 (cid:1)(cid:19) , n ≥ . (A7)Note that in some presentations, including [16], the first Chebyshev coefficient (cid:98) y is definedwith a 1/2 in front of it, so that the exceptions in (A7) for (cid:98) y and (cid:98) y can be written aspart of the general case. In either case, the formulation gives two free parameters (cid:98) y and (cid:98) y (cid:48) , which are obtained using the boundary conditions. To match the number of unknowns(coefficients) and equations, we assume (cid:98) y (cid:48)(cid:48) n = (cid:98) y (cid:48) n = 0 for n > N − (cid:98) y n using(A7).We can reformulate the boundary value problem using the Chebyshev series representa-6tions (A3) as N − (cid:88) n =0 ( (cid:98) y (cid:48)(cid:48) n − k (cid:98) y n ) T n ( z ) = N − (cid:88) n =0 (cid:98) f n T n ( z ) . (A8)Matching modes gives a system of equations for the Chebyshev coefficients (cid:98) y (cid:48)(cid:48) n − k (cid:98) y n = (cid:98) f n n = 0 , . . . N − , (A9) N − (cid:88) n =0 ( (cid:98) y (cid:48) n + k (cid:98) y n ) = α, N − (cid:88) n =0 ( (cid:98) y (cid:48) n − k (cid:98) y n )( − n = β. (A10)with N +2 equations and N +2 unknowns (cid:98) y , (cid:98) y (cid:48) , (cid:98) y (cid:48)(cid:48) , . . . (cid:98) y (cid:48)(cid:48) N − . We solve the algebraic system ofequations (A9) and (A10) for the second derivative coefficients (cid:98) y (cid:48)(cid:48) , . . . (cid:98) y (cid:48)(cid:48) N − and integrationconstants (cid:98) y and (cid:98) y (cid:48) , then determine (cid:98) y n by applying the double integral operation on thecoefficients (A7).If the domain is z ∈ [0 , H ], then the equations are (cid:98) y (cid:48)(cid:48) n − k H (cid:98) y n = (cid:98) f n n = 0 , . . . N − , (A11) N − (cid:88) n =0 (cid:18) H (cid:98) y (cid:48) n + k H (cid:98) y n (cid:19) = α, N − (cid:88) n =0 (cid:18) H (cid:98) y (cid:48) n − k H (cid:98) y n (cid:19) ( − n = β. (A12)For the k = 0 mode, the system reduces to N − (cid:88) n =0 (cid:98) y (cid:48)(cid:48) n T n (cid:18) zH − (cid:19) = N − (cid:88) n =0 (cid:98) f n T n (cid:18) zH − (cid:19) , (A13)which gives the trivial system of equations (cid:98) y (cid:48)(cid:48) n = (cid:98) f n n = 0 , . . . N − . (A14)Our solver uses homogeneous boundary conditions for the k = 0 mode, N − (cid:88) n =0 (cid:98) y n = 0 , N − (cid:88) n =0 (cid:98) y n ( − n = 0 . (A15)We use a Schur complement approach to solve the algebraic system of equations (A11)7and (A12). We can write the system in block form as A BC D (cid:98) y (cid:48)(cid:48) (cid:98) y (cid:98) y (cid:48) = (cid:98) f αβ . (A16)Here B is N × C is 2 × N , D is 2 ×
2, and A is an N × N pentadiagonal matrixwith only three nonzero diagonals, which we first (pre)factorize using a fast algorithm forsuch matrices ([39], specialized to the case of only three nonzero diagonals). Our Schurcomplement approach is then to solve the 2 × CA − B − D ) (cid:98) y (cid:98) y (cid:48) = CA − f − αβ (A17)for (cid:98) y and (cid:98) y (cid:48) . We then obtain the coefficients (cid:98) y (cid:48)(cid:48) = ( (cid:98) y (cid:48)(cid:48) , . . . (cid:98) y (cid:48)(cid:48) N − ) from back-substitution (cid:98) y (cid:48)(cid:48) = A − f − B (cid:98) y (cid:98) y (cid:48) . (A18) [1] Leslie Greengard. Fast algorithms for classical physics. Science , 265(5174):909–914, 1994.[2] DS Shamshirgar, R Yokota, A-K Tornberg, and Berk Hess. Regularizing the fast multipolemethod for use in molecular simulation.
The Journal of Chemical Physics , 151(23):234113,2019.[3] Jiuyang Liang, Jiaxing Yuan, Erik Luijten, and Zhenli Xu. Harmonic surface mapping algo-rithm for molecular dynamics simulations of particle systems with planar dielectric interfaces.
The Journal of Chemical Physics , 152(13):134109, 2020.[4] S Alireza Ghasemi, Alexey Neelov, and Stefan Goedecker. A particle-particle, particle-densityalgorithm for the calculation of electrostatic interactions of particles with slablike geometry.
J. Chem. Phys. , 127(22):224102, 2007.[5] Franziska Nestler, Michael Pippig, and Daniel Potts. Fast ewald summation based on nfftwith mixed periodicity.
Journal of Computational Physics , 285:280–315, 2015. [6] Dag Lindbo and Anna-Karin Tornberg. Spectral accuracy in fast ewald-based methods forparticle simulations. Journal of Computational Physics , 230(24):8744–8761, 2011.[7] Davood Saffar Shamshirgar and Anna-Karin Tornberg. Fast ewald summation for electrostaticpotentials with arbitrary periodicity. arXiv preprint arXiv:1712.04732 , 2017.[8] D. R. Ladiges, S. P. Carney, A. Nonaka, K. Klymko, G. Moore, A. L. Garcia, S. R. Natesh,A. Donev, , and J. B. Bell. A Discrete Ion Stochastic Continuum Overdamped Solvent Algo-rithm for Modeling Electrolytes. Submitted to Phys. Rev. Fluids, ArXiv preprint 2007.03036,2020.[9] J.R. Phillips and J.K. White. A precorrected-FFT method for electrostatic analysis of com-plicated 3-D structures.
IEEE Trans. Computer-Aided Design , 16(10):1059–1072, 1997.[10] Sandeep Tyagi, Axel Arnold, and Christian Holm. Electrostatic layer correction with imagecharges: A linear scaling method to treat slab 2 d+ h systems with dielectric interfaces.
TheJournal of chemical physics , 129(20):11B616, 2008.[11] Alexandre P dos Santos, Matheus Girotto, and Yan Levin. Simulations of coulomb systemsconfined by polarizable surfaces using periodic green functions.
The Journal of chemicalphysics , 147(18):184105, 2017.[12] Alexandre P dos Santos and Yan Levin. Electrolytes between dielectric charged surfaces:Simulations and theory.
The Journal of Chemical Physics , 142(19):194104, 2015.[13] Axel Arnold, Konrad Breitsprecher, Florian Fahrenberger, Stefan Kesselheim, Olaf Lenz, andChristian Holm. Efficient algorithms for electrostatic interactions including dielectric con-trasts.
Entropy , 15(11):4569–4588, 2013.[14] Dag Lindbo and Anna-Karin Tornberg. Fast and spectrally accurate ewald summation for2-periodic electrostatic systems.
The Journal of chemical physics , 136(16):164111, 2012.[15] Anna-Karin Tornberg. The ewald sums for singly, doubly and triply periodic electrostaticsystems.
Advances in Computational Mathematics , 42(1):227–248, 2016.[16] Leslie Greengard. Spectral integration and two-point boundary value problems.
SIAM Journalon Numerical Analysis , 28(4):1071–1080, 1991.[17] L. Greengard and J. Lee. Accelerating the nonuniform fast fourier transform.
SIAM Review ,46(3):443–454, 2004.[18] Alexander H Barnett, Jeremy Magland, and Ludvig af Klinteberg. A parallel nonuniform fastfourier transform library based on an “exponential of semicircle” kernel.
SIAM Journal on Scientific Computing , 41(5):C479–C504, 2019.[19] Lloyd N Trefethen.
Spectral methods in MATLAB , volume 10. Siam, 2000.[20] Mark D Tinkle and SE Barlow. Image charge forces inside conducting boundaries.
Journal ofApplied Physics , 90(3):1612–1624, 2001.[21] Nawaf Bou-Rabee and Eric Vanden-Eijnden. Pathwise accuracy and ergodicity of metropolizedintegrators for sdes.
Communications on Pure and Applied Mathematics , 63(5):655–696, 2010.[22] N. Bou-Rabee, A. Donev, and E. Vanden-Eijnden. Metropolis Integration Schemes for Self-Adjoint Diffusions.
SIAM J. Multiscale Modeling and Simulation , 12(2):781–831, 2014.[23] B Leimkuhler, C Matthews, and MV Tretyakov. On the long-time integration of stochasticgradient systems. In
Proc. R. Soc. A , volume 470, page 20140120. The Royal Society, 2014.[24] V Ballenegger, A Arnold, and JJ Cerda. Simulations of non-neutral slab systems with long-range electrostatic interactions in two-dimensional periodic boundary conditions.
The Journalof chemical physics , 131(9):094107, 2009.[25] Raul P. Pelaez. Uammd. https://github.com/RaulPPelaez/UAMMD , 2020.[26] Benedict Leimkuhler, Charles Matthews, and Gabriel Stoltz. The computation of averagesfrom equilibrium and nonequilibrium Langevin molecular dynamics.
IMA Journal of Numer-ical Analysis , 36(1):13–79, 01 2015.[27] T. Croxton, D. A. McQuarrie, G. N. Patey, G. M. Torrie, and J. P. Valleau. Ionic solution nearan uncharged surface with image forces.
Canadian Journal of Chemistry-Revue CanadienneDe Chimie , 59(13):1998–2003, 1981.[28] Zhenli Xu, Manman Ma, and Pei Liu. Self-energy-modified poisson-nernst-planck equations:Wkb approximation and finite-difference approaches.
Physical Review E , 90(1):013307, 2014.[29] Manman Ma, Zhenli Xu, and Liwei Zhang. Modified poisson-nernst-planck model withcoulomb and hard-sphere correlations. arXiv preprint arXiv:2002.07489 , 2020.[30] David Andelman. Electrostatic properties of membranes: the poisson-boltzmann theory. In
Handbook of biological physics , volume 1, pages 603–642. Elsevier, 1995.[31] Amir Levy, J Pedro de Souza, and Martin Z Bazant. Breakdown of electroneutrality innanopores.
Journal of Colloid and Interface Science , 2020.[32] A. M. Fiore, F. Balboa Usabiaga, A. Donev, and J. W. Swan. Rapid sampling of stochas-tic displacements in brownian dynamics simulations.
J. Chem. Phys. , 146(12):124116, 2017.Software available at https://github.com/stochasticHydroTools/PSE . [33] Dag Lindbo and Anna-Karin Tornberg. Spectrally accurate fast summation for periodic stokespotentials. Journal of Computational Physics , 229(23):8994–9010, 2010.[34] J. P. Hernandez-Ortiz, J. J. de Pablo, and M. D. Graham. Fast Computation of Many-Particle Hydrodynamic and Electrostatic Interactions in a Confined Geometry.
Phys. Rev.Lett. , 98(14):140602, 2007.[35] Shriram Srinivasan and Anna-Karin Tornberg. Fast ewald summation for green’s functions ofstokes flow in a half-space.
Research in the Mathematical Sciences , 5(3):35, 2018.[36] Wen Yan and Michael Shelley. Universal image systems for non-periodic and periodic stokesflows above a no-slip wall.
Journal of Computational Physics , 375:263–270, 2018.[37] Wen Yan and Robert Blackwell. Kernel aggregated fast multipole method: Efficient summa-tion of laplace and stokes kernel functions. arXiv preprint arXiv:2010.15155 , 2020.[38] Z. Gimbutas, L. Greengard, and S. Veerapaneni. Simple and efficient representations for thefundamental solutions of Stokes flow in a half-space.
Journal of Fluid Mechanics , 776:R1,2015. Code available at .[39] AA Karawia. Two algorithms for solving general backward pentadiagonal linear systems.