A Discrete Ion Stochastic Continuum Overdamped Solvent Algorithm for Modeling Electrolytes
Daniel R. Ladiges, Sean P. Carney, Andrew Nonaka, Katherine Klymko, Guy C. Moore, Alejandro L. Garcia, Sachin R. Natesh, Aleksandar Donev, John B. Bell
AA Discrete Ion Stochastic Continuum Overdamped Solvent Algorithm for ModelingElectrolytes
D. R. Ladiges, ∗ A. Nonaka, K. Klymko, G. C. Moore, and J. B. Bell
Center for Computational Sciences and Engineering, LBNL † S. P. Carney
Department of Mathematics, UT Austin
A. L. Garcia
Department of Physics and Astronomy, SJSU
S. R. Natesh and A. Donev
Courant Institute of Mathematical Sciences, New York University, New York, NY 10012 (Dated: August 25, 2020)In this paper we develop a methodology for the mesoscale simulation of strong electrolytes. Themethodology is an extension of the Fluctuating Immersed Boundary (FIB) approach that treats asolute as discrete Lagrangian particles that interact with Eulerian hydrodynamic and electrostaticfields. In both algorithms the Immersed Boundary (IB) method of Peskin is used for particle-fieldcoupling. Hydrodynamic interactions are taken to be overdamped, with thermal noise incorporatedusing the fluctuating Stokes equation, including a “dry diffusion” Brownian motion to account forscales not resolved by the coarse-grained model of the solvent. Long range electrostatic interactionsare computed by solving the Poisson equation, with short range corrections included using a novelimmersed-boundary variant of the classical Particle-Particle Particle-Mesh (P3M) technique. Alsoincluded is a short range repulsive force based on the Weeks-Chandler-Andersen (WCA) potential.The new methodology is validated by comparison to Debye-H¨uckel theory for ion-ion pair correlationfunctions, and Debye-H¨uckel-Onsager theory for conductivity, including the Wien effect for strongelectric fields. In each case good agreement is observed, provided that hydrodynamic interactionsat the typical ion-ion separation are resolved by the fluid grid. ∗ [email protected] † https://ccse.lbl.gov a r X i v : . [ phy s i c s . c o m p - ph ] A ug I. INTRODUCTION
Understanding transport properties in electrolytes is important for the study of fundamental processes such aselectrophoresis, electro-osmosis and electrochemistry that arise in both biological systems [1] and engineered devices[2] such as catalytic micropumps [3, 4], batteries [5, 6], and fuel cells [7–9]. Many of these phenomena occur at themesoscale, where the effects of thermal fluctuations must be captured correctly. This has previously motivated thedevelopment of fluctuating hydrodynamics [10] methods for electrolytic flows [11–13], in which the system is describedusing stochastic partial differential equations, which are solved on a grid. Using this approach, thermal fluctuationsare captured while avoiding the computational expense inherent in direct molecular simulations [14]. However, onedrawback of a purely continuum model is that it is not suitable for mesoscopic systems that contain features that occuron a molecular scale, but which still have an important effect on the overall dynamics. In the context of electrolytesthe electrical double layer [15] effect is an important example of such phenomena. In order to capture this behavior while avoiding the computational expense of molecular dynamics, an approachis required that combines a particle model for the solutes with a coarse-grained model of the solvent. Currently, themost well established numerical method that uses this approach is Brownian dynamics (BD) [16]. The key feature ofthis technique is the use of Green’s functions to represent particles immersed in an implicit solvent. These Green’sfunctions describe the effects of single particles, which can be joined together in groups to form complex moleculessuch as polymers, or even structures such as flexible membranes. While this approach works well in many cases, ithas a range of shortcomings. First, in a basic implementation, the calculation of inter-particle hydrodynamic forceshas computational complexity O ( N ), where N is the number of particles being simulated. Further, the calculationof the thermal noise applied to each particle scales as O ( N ) (from the Cholesky factorization of the mobility matrix),making simulation of large systems infeasible. Additionally, the application of boundary conditions is complicated toimplement in all but the simplest cases.More recently, an approach referred to as General Geometry Ewald-like Method (GGEM) [17, 18] has been proposed.Similar in principle to the Particle-Particle Particle-Mesh (P3M) method [19] used in electrostatics, short rangehydrodynamic interactions are computed using Green’s functions, and long range interactions are computed usinga grid-based discretization of the relevant PDE, in this case Stokes’ equation. With suitable preconditioning of theStokes solver and choice of the splitting parameter, this achieves near linear scaling for deterministic hydrodynamics bylimiting the calculation of pairwise interactions to small subsets of particles. In principle, GGEM allows for handling ofnontrivial boundary conditions, although some approximations are introduced [17, 18] to handle the boundaries in thenear field; see Ref. 20 for exact handling of a bottom wall geometry with Ewald splitting. Importantly, the handlingof Brownian motion in GGEM is based on methods introduced in the 1970s by Fixman, and these methods increasethe computational cost over deterministic simulations many-fold [21]. First, generating the Brownian displacementsuses an iterative method, which ideally employs preconditioners [22, 23] in order to control the number of iterations.Second, when boundaries are present, generating the correct stochastic drift required to achieve discrete fluctuation-dissipation balance relies on a method of Fixman that requires generating the action of the square root of the resistancematrix, which itself requires a slowly-converging iterative method.A related approach that addresses these issues is the Fluctuating Immersed Boundary (FIB) [24] method. Thistechnique uses Immersed Boundary (IB) [25] kernels instead of Green’s functions to couple particles to an explic-itly simulated solvent. Importantly, whereas both the original IB and GGEM methods use standard deterministichydrodynamics to represent the solvent, in FIB the solvent is simulated using an incompressible fluctuating hydro-dynamics methodology [26]. Using this approach, thermal noise arises directly from the fluid, avoiding the complexper particle calculations required by GGEM and other Brownian dynamics methods. In FIB, no iterative methodsare required to generate either the Brownian increments or the stochastic drift (rather, they are computed togetherwith the deterministic displacements or generated by a simple “random finite difference” midpoint scheme inspiredby but different from Fixman’s), and discrete fluctuation-dissipation balance is ensured exactly by construction evenin the presence of boundaries. Linear time complexity, in both domain size and particle number, is achieved for thecalculation of hydrodynamic interactions and Brownian motion under all conditions. Additionally, the application ofconfining boundary conditions, such as no-slip walls, is trivial when using FIB. We note a similar approach is taken inthe Stochastic Eulerian Lagrangian Method (SELM), described in Ref. 27; some comments regarding the relationshipbetween FIB and SELM are given in Ref. 24.Recently, the FIB approach to Brownian motion has been combined with the Ewald splitting used in GGEM forperiodic systems in the Positively Split Ewald (PSE) method [21]; we note that the Hasimoto splitting used in PSE is For example, the thickness of the double layer is roughly one Debye length, λ D , and for a molar concentration of 1.0M there is only oneion per (3 λ D ) . In particular, in GGEM, the boundary condition imposed in the Stokes solver involves the near field, which varies on scales that cannot be resolved by the grid. A further point-particle approximation is made in the boundary conditions to keep the mobility positivesemi-definite [17, 18]. the Fourier space equivalent of the real-space splitting used in GGEM [28]. PSE is a linear-scaling, spectrally-accurate,and grid-independent method that has been shown [21] to be an order of magnitude faster than traditional approachesbased on the techniques developed by Fixman. However, including boundaries in spectral-Ewald methods relying onFourier transforms such as PSE is nontrivial and has not, to our knowledge, been done yet even for deterministicsimulations.In this paper, we extend the FIB method to the simulation of electrolytic flows. Using this approach, individualcharged ions are described using immersed boundary kernels, and the solvent is treated using fluctuating hydrody-namics. In doing so we incorporate two new contributions to the FIB methodology: (i) a “dry diffusion” process whichimplements a coarse-grained model for small ions (e.g., Na + ), based on the approach derived in Ref. 29, and (ii) theincorporation of electrostatic effects with an immersed-boundary P3M implementation for efficiently and accuratelycomputing the electric field. We refer to this approach as Discrete Ion Stochastic Continuum Overdamped Solvent(DISCOS).The layout of this paper is as follows. In Sec. II, we summarize the Brownian dynamics approach and describe itsrelationship to the FIB method. In Sec. II E we discuss the dry diffusion process. Sec. III contains a description ofelectrostatic and close range forces, while Sec. IV details the numerical methodology used to implement the DISCOSalgorithm. In Sec. V the DISCOS method is tested by comparison to theoretical results for the radial distributionfunction and electrical conductivity. This includes an analysis of the effect of the dry diffusion approach described inSec. II E. Finally, some concluding remarks are given in Sec. VI. Note that in this article, to fully detail and test thenovel aspects of the method without introducing additional complications, we describe the application of DISCOS tounbounded/periodic systems. The inclusion of boundaries will be the subject of a future publication. II. STOCHASTIC HYDRODYNAMICS
This section formulates stochastic hydrodynamics for Brownian particles, and Sec. III describes the additionalintermolecular forces present when these particles are ions. First we briefly review the Brownian dynamics andfluctuating immersed boundary approaches to simulating particle-fluid systems; more detailed descriptions can befound in Refs. 16, 17, 24, 30 and 31. Both approaches assume an infinite Schmidt number,Sc = ηρD → ∞ , (1)where D is the diffusion coefficient of a particle, and η and ρ are the viscosity and density of the fluid, respectively. Inthis asymptotic regime, the diffusion time of the particles is large compared to the relaxation time of the fluid, and theflow can therefore be treated as quasi-steady. A detailed discussion of the relationship between particle diffusion andSchmidt number is given in Ref. 32, and the validity of this approximation in the context of electrolytes is discussedin Sec. V.The FIB method was originally intended as a fast algorithm for performing Brownian dynamics with hydrodynamicinteractions for colloidal suspensions; the reader can consult Refs. 21, 33 for related state-of-the-art spectrally-accuratemethods for periodic suspensions. Here we adapt this methodology to ionic solutions. This necessitates some importantchanges to account for the fact that ions are much smaller than typical colloids, as we describe in detail in SectionsII B and II E. A. Brownian dynamics with hydrodynamic interactions
We consider a solute (ions) represented as N particles with positions x ( t ) = { x , · · · , x i , · · · , x N } . The equationof motion for the particles is d x dt = M F + (cid:112) k B T M / (cid:5) W (2)= M F + (cid:112) k B T M ◦ M − / W . (3)Here, F ( x , t ) = { F , · · · , F i , · · · , F N } are the net forces acting on the particles (in this article they consist of short-ranged intermolecular as well as long-ranged electrostatic forces), W ( t ) = { W , · · · , W i , · · · , W M } are independentGaussian white noise processes, k B is Boltzmann’s constant, and T is the temperature. The stochastic product symbol Note that it is not required that M = N . (cid:5) indicates the kinetic interpretation [34, 35] of the stochastic integral, which can be viewed a mixed Stratonovich-Itˆointerpretation, with ◦ denoting the Stratonovich and the absence of a symbol denoting the Itˆo stochastic product. Thesymmetric positive-definite mobility matrix M ( x ) encodes the hydrodynamic interactions between particles. Here,we define M / and M − / to be (not necessarily square) matrices such that M = M / (cid:16) M / (cid:17) (cid:63) , M − = M − / (cid:16) M − / (cid:17) (cid:63) , (4)where the (cid:63) indicates L adjoint (conjugate-transpose). When the kinetic interpretation of the stochastic integral isapplied this property ensures fluctuation-dissipation balance is obeyed, i.e., that the equilibrium dynamics are timereversible with respect to the Gibbs–Boltzmann distribution. When we recast Eq. (2) as a system of Itˆo stochasticdifferential equations for numerical simulation, it is necessary to include a stochastic drift term, yielding d x dt = M F + k B T ∇ x · M + (cid:112) k B T M / W , (5)where ∇ x · is the divergence with respect to the particle position variables.In Brownian dynamics the main considerations in the numerical integration of Eq. (5) are the construction and fastapplication of M , and the subsequent application of M / and calculation of ∇ x · M . The action of a point forceat the origin on the fluid is given by the Green’s function solution to the steady Stokes’ equation, the Oseen tensor O ( r ) = 18 πηr (cid:18) I + r ⊗ r r (cid:19) . (6)Here, r is the position vector with r = | r | and I is the identity tensor. This can be modified to include finite size [36]and close range corrections, giving the Rotne-Prager-Yamakawa (RPY) tensor [37], R ( r ; a ) = 16 πηa C I + C r ⊗ r r r > aC I + C r ⊗ r r r ≤ a , (7)with C = 3 a r + a r , C = 3 a r − a r ,C = 1 − r a , C = 3 r a , (8)where we have assumed all particles have the same hydrodynamic radius a . The 3 × i due to a force applied to particle j is then M ij = R ( x ij ; a ) , (9)where x ij = x i − x j . The complete mobility matrix is assembled by computing Eq. (9) for all particle pairs. We notethat M is symmetric positive-definite by construction. Having obtained M , the matrix M / can be calculatedexactly using Cholesky decomposition or applied using a more rapid iterative method of approximation [22, 38].Rather than directly calculating the divergence of the mobility matrix to obtain the stochastic drift, the driftis traditionally included via the Fixman [16, 39] midpoint time stepping algorithm. This can also be viewed as adirect discretization of the kinetic version of the stochastic integral [35] or, equivalently, a direct discretization of themixed-integral in Eq. (3): x n +1 / ,(cid:63) = x n + ∆ t M ( x n ) F n + (cid:114) ∆ tk B T M n ) / W n , x n +1 = x n + ∆ t M ( x n +1 / ,(cid:63) ) (cid:32) F n +1 / ,(cid:63) + (cid:114) k B T ∆ t (cid:16) M ( n ) (cid:17) − / W n (cid:33) . (10)Here, ∆ t is the discrete time step and W is a vector of independent Gaussian random variables. The superscripts n and n + 1 denote values at consecutive time steps, with n + 1 / , (cid:63) indicating a midpoint update. B. Fluctuating Hydrodynamics Formulation
As discussed in Sec. I, the FIB method differs from Brownian dynamics principally in that it explicitly simulates acoarse grained model of the solvent [29] on an Eulerian grid. The solvent is taken to be isothermal and incompressible,and is therefore modeled by the fluctuating Stokes equations [29, 40], ∂ v ∂t + ∇ r p − η ∇ r v = f + (cid:112) k B T η ∇ r · Z , (11a) ∇ r · v = 0 , (11b)where ∇ r is the gradient operator with respect to r , v ( r , t ) is the fluid velocity, p ( r , t ) is the pressure, and f ( r , t ) is aforce density applied to the fluid; this is the mechanism by which the immersed Brownian particles interact with thesolvent. Finally, Z ( r , t ) is a random, symmetric, spatio-temporal Gaussian tensor field whose components are whitein space and time, i.e., they have mean zero and covariances (cid:104)Z ij ( r , t ) Z kl ( r (cid:48) , t (cid:48) ) (cid:105) = ( δ ik δ jl + δ il δ jk ) δ ( t − t (cid:48) ) δ ( r − r (cid:48) ) (12) ≡ Σ δ ( t − t (cid:48) ) δ ( r − r (cid:48) ) . Here, δ ij is the Kronecker delta function, δ is the Dirac delta function, Σ is the pointwise covariance of the noise,and (cid:104)· · · (cid:105) denotes an ensemble average. This time-dependent system is a set of linear, stochastic partial differentialequations with additive noise, and it therefore has a well-defined mathematical interpretation [41]. Furthermore, itcan be shown to satisfy fluctuation-dissipation balance [27, 29].In Refs. 24, 27, 40 the diffusing particles are coupled directly to the continuum equations, Eq. (11). However,as discussed in detail in Refs. 29, 42, fluctuating hydrodynamics (FHD) is a coarse-grained description that onlyhas a clear interpretation once the conserved fluid quantities (mass, momentum, energy) are spatially coarse-grainedon a discrete grid. Each hydrodynamic cell should contain a sufficiently large number of solvent molecules for thecoarse-grained description to be justified. Since our particles are (solvated) ions of size comparable to the solventmolecules, this implies that the fluid grid size must be larger than a diffusing “nano-particle”. This is exactly the caseanalyzed in Ref. 29 using the theory of coarse-graining, where the discrete FHD equations are derived by coupling asingle nano-particle with compressible, isothermal FHD. We base our discrete equations on this derivation, and takeadvantage of the flexibility of the framework to improve numerical properties, such as grid-independence of physicalobservables and preservation of translational invariance.One key conclusion of the analysis in Ref. 29 is that the total or effective diffusion coefficient of an immersed particleconsists of two pieces. The first comes from the random advection of the particle by the coarse-grained stochasticfluid velocity; we will refer to this process as “wet diffusion”, with the diffusion coefficient D wet . The second comesfrom the additional random motion of the particle relative to the coarse-grained fluid velocity; we will term this“dry diffusion”, with the coefficient D dry . One can formally write a Green-Kubo formula expressing the dry diffusioncoefficient as the integral of the autocorrelation function of the particle velocity relative to the coarse-grained fluidvelocity [29]; we discuss how we set the value of the dry diffusion coefficient in Sec. II E. C. Steady, overdamped limit
For high Schmidt numbers the fluid relaxation time is fast on the time scale of particle diffusion, so the time-dependent system Eq. (11) can be accurately approximated by a (quasi-)steady Stokes system. This limiting behavior,referred to as the overdamped limit, can be derived by a rigorous analysis of the limit as Sc → ∞ ; this is shown inAppendix A of Ref. 40. In the overdamped limit, the dynamic variables are just the positions of the immersedparticles, and they follow an equation that is identical in structure to the overdamped Langevin equation, Eq. (2).Taking the subscript h to indicate a discrete operator, in the overdamped limit we take the discrete fluid velocityand pressure to satisfy the steady-state form of Eq. (11) given by ∇ h p − η ∇ h v = f + (cid:114) k B T η ∆ V ∇ h · Z , (13a) ∇ h · v = 0 , (13b)where Z ( r h , t ) is a finite dimensional collection of white noise processes representing the spatial discretization of Z on a regular grid with positions r h , and ∆ V is the cell volume. Preserving the fluctuation-dissipation property thenrequires that the discrete divergence be the negative of the L adjoint of the discrete gradient, ∇ (cid:63)h = − ∇ h · , and that ∇ h = ∇ h · Σ ∇ h = − ∇ (cid:63)h Σ ∇ h . (14)Note that the spatial discretization scheme described in Sec. IV below is designed to ensure that these properties hold. D. Particle-fluid coupling
In the FIB method the spatial extent of each particle is defined by a compact kernel function δ hy ( r ), where thesuperscript “hy” indicates that this kernel applies to hydrodynamic fields (the analogous electrostatic kernel will bedenoted δ es ). The kernel is used to interpolate the local fluid velocity to the particles’ locations x i , and to definethe region over which force on the particle is transmitted (spread) to the fluid; the functional form of the kernelis discussed in Section IV A. To perform these operations, we define discrete interpolation and spreading operators, J hy h ( x i ) and S hy h ( x i ), such that V i = J hy h ( x i ) v = (cid:90) δ hy ( x i − r h ) v ( r h ) d r h , (15) f s ( r h ) = (cid:16) S hy h ( x i ) F (cid:17) ( r h ) = N (cid:88) i =1 δ hy ( x i − r h ) F i , (16)where the hydrodynamic (advective) part of the particle velocity, i.e., the velocity of the fluid at the location of theparticles, is denoted by V = { V , · · · , V i , · · · , V N } . The force density spread to the fluid from the particles is givenby f s , and as above, we use the subscript h to indicate discrete operators. Note that the integral in Eq. (15) is justa short-hand notation for a weighted sum over the grid points; the exact form of the discretization is discussed inSec. IV A.It is important to note that this coupling of the particles to the coarse-grained fluid, taken from immersed boundarymethodology, is different from (though closely related to) that used in the theory of coarse-graining given in Ref. 29.The FIB formulation is advantageous because it is simpler to implement, achieves much better translational invariancein practical numerical implementations, and allows flexibility in choosing δ hy .In addition to advection by the velocity given by Eq. (15), particles move by the dry diffusion process, giving thecomplete overdamped equations of motion for the DISCOS method (see also Appendix B.1 in Ref. 43), d x i dt = V i + k B T ( ∇ x · M ) i wet + D dry i k B T F i + (cid:113) D dry i W dry i dry , (17)where W dry i ( t ) is another independent Gaussian white noise process. Note that here we have written the equation ofmotion for a single particle, since D dry i depends on the particle’s species.We define the discrete Stokes operator L h such that the solution of Eq. (13) is v = L − h (cid:32) f + (cid:114) k B T η ∆ V ∇ h · Z (cid:33) , (18)so that the particles’ hydrodynamic velocity is V i = J hy h ( x i ) L − h (cid:32) S hy h ( x i ) F ( x ) + (cid:114) k B T η ∆ V ∇ h · Z (cid:33) . (19)This allows us to relate the DISCOS method to Brownian dynamics (Eq. (2)) by observing that M = Diag (cid:40) D dry k B T (cid:41) + J hy h L − h S hy h , (20)where D dry = (cid:110) D dry1 , · · · , D dry i , · · · , D dry N (cid:111) , and Diag { X } denotes a diagonal matrix with the values of X on thediagonal. The Brownian velocity is expressed as M / W ≡ Diag (cid:115) D dry k B T W dry i + (cid:114) η ∆ V J hy h L − h ∇ h · Z , (21)where (cid:112) D dry = (cid:26)(cid:113) D dry1 , · · · , (cid:113) D dry i , · · · , (cid:113) D dry N (cid:27) . In Ref. 24 (see also Appendix B of Ref. 43) it is demonstratedthat this combination satisfies fluctuation-dissipation balance, provided Eq. (14) holds, and the spreading and inter-polation are related by the adjoint property (cid:16) J hy h (cid:17) (cid:63) = ∆ V S hy h , which is assured by the immersed-boundary method.What is now needed is a way to evaluate the stochastic drift. Direct implementation of the Fixman method(Eq. (10)) is complicated by the fact that it requires the inverse of the mobility matrix, which is difficult to obtainin the context of FIB. As demonstrated in Ref. 24, the divergence of the mobility can be evaluated numerically usinga finite difference method. However this would require multiple evaluations of M per step, which requires multipleapplications of L − h , i.e., multiple solves of Stokes equation. This can be avoided by re-writing the stochastic driftusing the chain rule ∇ x · M = J hy h L − h ∇ x · S hy h + (cid:16) ∇ x J hy h (cid:17) : (cid:16) L − S hy h (cid:17) (22)where, in summation notation, the double dot product is defined as (cid:110)(cid:16) ∇ x J hy h (cid:17) : (cid:16) L − h S hy h (cid:17)(cid:111) ij = ( ∂ j ( J hy h ) ik )( L h ) − kl ( S hy h ) lj . (23)Defining the thermal forcing f th = k B T ∇ x · S hy h , (24)the first term on the right hand side of Eq. (22) is then accounted for by including f th as part of the force density inEq. (13), f = f s + f th . (25)The divergence of the spreading operator can be evaluated with a single application of L − h ; the details are discussedin Sec. IV. The second part of the stochastic drift (corresponding to the second term on the right hand side of Eq. (22))is obtained by a midpoint temporal integrator described in Sec. IV D, similar in spirit to how the Fixman midpointscheme (Eq. (10)) obtains the total stochastic drift. E. Total, wet, and dry diffusion
At the end of Sec. II B we describe two diffusion processes: wet diffusion, which arises from fluctuations in thecoarse grained hydrodynamic solution, and dry diffusion, which represents additional Brownian motion that is notresolved by the coarse graining process. Here we explain how we set the dry diffusion coefficient.Consider a single freely diffusing isolated spherical nanoparticle suspended in a quiescent fluid. The particle willperform a standard Brownian motion with a “total” diffusion coefficient D tot . From the Stokes-Einstein relation, thetotal diffusion coefficient for a sphere with radius a t is D tot = k B Tςηa t , (26)where ς is a constant depending on the boundary condition on the particle; here we take ς = 6 π , corresponding to ano slip boundary condition. Note that as long as D tot is held fixed the relative value of ς and a t will have no bearingon the dynamics.As described in Sec. II B, when simulating such a particle using the FIB method the fluctuations in the hydrodynamicvelocity will yield a wet diffusion process with coefficient D wet . Because these fluctuations are transmitted to theparticle via the interpolation operator, J hy h , the effective hydrodynamic radius of the particle will depend on the formof the kernel function δ hy ; we designate this wet radius a w , where D wet = k B Tςηa w . (27)The value of a w is determined by the grid used to solve Eq. (13) – coarsening this grid will result in an increase of a w and a reduction in the diffusion experienced by the particle.If the discretization is such that a w > a t , then the additional dry diffusion process to be added is such that D dry = k B Tςηa d , (28)where the dry radius a d is a d = a w a t a w − a t , (29)giving D tot = D wet + D dry . (30)In order for the dry diffusion to remain positive, the total radius a t must be smaller than the wet radius a w . Inpractice, D wet is set by the mesh spacing and the choice of kernel. We then specify D dry so that each species diffuseswith the desired total diffusion coefficient. Using this approach, the diagonal terms in an under-resolved mobilitymatrix are corrected with the dry diffusion terms; see Eq. (20). This has the effect of greatly improving computationalspeed, at the cost of neglecting short range hydrodynamic contributions. For electrolytes, in many cases we expectthese contributions to be small compared to the electrostatic effects; this is discussed further in Section V B.The above procedure was suggested, though not implemented, in Ref. 29, where the total, wet, and dry diffusioncoefficients are referred to as the renormalized, enhancement, and bare diffusion coefficients, respectively. A similarapproach is also used for inertial systems in Refs. 43–45. Note that here we have taken the dry diffusion coefficient tobe isotropic and spatially homogeneous, which is true only in the case of unbounded/periodic domains. In this case,no stochastic drift correction is necessary for the dry component of Eq. (17). For other boundary conditions the drydiffusion coefficients can be spatially dependent, and a corresponding stochastic drift term would have to be included. III. IONIC FORCES
As discussed in Sec. I, we are considering long range electrostatic and short range repulsive forces acting on theions. For a given ion i , the total force is given by F i = F E i + (cid:88) j ∈ Ω R i F R ij + F ext i , (31)where F E i is the electrostatic force, F ext i indicates forces due to an applied field (e.g., gravity, external electric field), F R ij is the short range repulsive force between particles i and j , and Ω R i indicates all particles within a given range ofthe i th particle. A. Electrostatic Forces
First, we describe our immersed-boundary variant of the classical Particle-Particle, Particle-Mesh (P3M) approach[14, 19] in detail, since we are not aware of prior work that implements an IB-P3M algorithm. We employ anapproach that is essentially linear in system size to evaluate electrostatic interactions, while recovering a point chargerepresentation of the ions which cannot be captured using purely grid based methods. For periodic systems, thecurrent state-of-the-art are variants of the Spectral Ewald method described in Ref. 46; a key advantage of our IBapproach is the ease of handling other types of boundary conditions that are crucial for modeling confined systems,as we will explore in future work. Some additional numerical results are presented in Appendix A.In the quasi-electrostatic approximation assumed here, the electrostatic force is found by solving Poisson’s equationfor the electrical potential φ , ∇ r φ = − (cid:37)/(cid:15), (32)where the charge density is (cid:37) ( r ) = N (cid:88) i =1 δ ( x i − r ) q i , (33)with q i being the charge of ion i , and (cid:15) the electrical permittivity of the solvent. The resulting electric field is thengiven by E = − ∇ r φ, (34)with the electrostatic force on ion i being F E i = q i E ( x i ); note that the singular self-induced electric field is definedto be zero.The P3M approach entails solving for the electrical potential on a grid, and then making local close range correctionsto the force in order to treat the ions as point charges. Note that for triply-periodic domains, the Spectral Ewaldmethod described in Ref. 46 is the most accurate technique. However, this method is based on Fourier transformsand is nontrivial to generalize to other boundary conditions and domain configurations.In the initial stage of the electrostatic force computation we map the particle charges to a mesh, discretize andsolve Eq. (32) on the mesh, and interpolate the resulting field back to the particle locations. Analogous to thehydrodynamics, the Lagrangian particles interact with the Eulerian electrostatic mesh via interpolation and spreadingoperators. E P i = J es h ( x i ) E = (cid:90) δ es ( x i − r h ) E ( r h ) d r h , (35) (cid:37) ( r h ) = S es h ( x i ) q = N (cid:88) i =1 δ es ( x i − r h ) q i , (36)where the integral in Eq. (35) is a short-hand notation for a weighted sum over grid points, just as for the hydrodynamicinterpolation operator. Here E P = (cid:110) E P1 , · · · , E P i , · · · , E P N (cid:111) is the smoothed value of E at the particle locations, q = { q , · · · , q i , · · · , q N } is the vector of particle charges, and δ es is the kernel used for the electrostatic grid; we notethat δ es can be chosen independently from δ hy . The numerical solution of Eq. (32) is discussed in Sec. IV C. Thesmoothed electrostatic force on particle i is then given by F P i = q i E P i . (37)A local correction is needed because a point charge representation of the ions would require an unreasonably finemesh when solving Eq. (32); for finite, non-zero ∆ r , the force calculated between two particles, i and j , on the meshbecomes less accurate as the distance between the particles becomes smaller. This effect becomes particularly severewhen the kernels overlap and the particles no longer appear as points, i.e., x ij < ψ , where ψ is on the order of thediameter of the kernel. The P3M method accounts for such situations by replacing the Poisson solution with a directcalculation of the Coulomb force, F C ij = 14 π(cid:15) q i q j x ij ˆ x ij , (38)where x ij = | x ij | , and the unit vector ˆ x ij = x ij /x ij . Specifically, the electrostatic force is computed as F E i = F P i + (cid:88) j ∈ Ω E i F LC ij , (39)where Ω E i indicates all particles within a range of ψ from particle i . Here the local or near-field correction is F LC ij = F C ij − F P ij , (40)where F P ij is the force computed by the immersed-boundary method if only the charges i and j were present. In thiswork the cutoff ψ ≤ L/
2, where L is the length of the periodic domain, so that for each pair of charges the localcorrection is only applied for the nearest periodic image of one of the charges.Following previous P3M and (Spectral) Ewald methods, the local correction is approximated as a radially-symmetricshort-ranged force F LC ij ≈ q i q j π(cid:15) (∆ r ) F LC (cid:16) x ij ∆ r (cid:17) ˆ x ij , (41)where the dimensionless scalar function F LC ( y ) is tabulated for 0 ≤ y ≤ ψ/ (∆ r ) in a pre-computation step, which isdone once and only once for a given kernel δ es . This table is computed using two charges in a domain which is large0enough that boundary effects are negligible. Because there is a small degree of translational variance, the solution foreach separation distance is calculated many times using random placements and orientations of the pair relative tothe grid, and the results averaged. We show numerical results for the near-field corrections in Appendix A.It is important to emphasize a key difference between the IB approach described here and that used in other P3M-style methods, including the spectral Ewald method [46]. In the latter, the smearing/smoothing of the delta functionsin Eq. (33) is done at the continuum level, using a smooth kernel function such as a Gaussian. For an unboundeddomain, the correction potential F LC ij can easily be computed analytically using Fourier-space integration. The smoothcontinuum density is then spread to the grid on which the Poisson equation is solved (using, e.g., a spectral method inthe Spectral Ewald method or a finite-difference approach as done here), and the solution is interpolated back on theparticles just as we do. The spreading and interpolation operations to and from the grid are done to spectral accuracyusing the non-uniform FFT in the Spectral Ewald method, while in classical P3M methods it is done using a separate“window” or “charge assignment function,” such as a power of the sinc function [19]; a plethora of methods exist withdifferent choices of how this step is done [14]. By contrast, in our approach there is only one kernel, δ es , that is usedboth to smear the charge and to communicate between particles and grid. In Appendix A, we study the accuracy ofour IB-P3M approach and isolate the contributions to the dominant error coming from the loss of translational androtational invariance from the Poisson solver and from the spreading and interpolation using a compactly-supportedkernel. B. Steric Repulsion
To compute the short range repulsive force F R ij , we model the finite size and excluded volume of the ions using amollified Weeks-Chandler-Andersen (WCA) interaction potential [47], which is a shifted and truncated Lennard-Jonespotential. This has been selected to approximate hard sphere repulsion whilst remaining numerically tractable. It isgiven by U sr (˜ x ; σ, ξ ) = ξ (cid:18) σ ˜ x − σ ˜ x (cid:19) (˜ x − ˜ x m ) − ξ (cid:18) σ ˜ x − σ ˜ x (cid:19) + ξ, ˜ x ≤ ˜ x m ξ (cid:18)(cid:16) σ ˜ x (cid:17) − (cid:16) σ ˜ x (cid:17) (cid:19) + ξ, ˜ x m < ˜ x < / σ , / σ ≤ ˜ x (42)where ξ is the magnitude of the potential, σ is the van der Waals diameter, and ˜ x is the radial distance between theparticles. Note that this potential has been selected to approximate a is set to produce a constant force for ˜ x ≤ ˜ x m equal to that for ˜ x = ˜ x m in order to reduce numerical stiffness; specific values of these parameters are discussed inSection V. The cutoff at 2 / σ ensures that only the repulsive part of the potential is used and shifting the potentialby ξ ensures that there is no discontinuity in the potential and in the force at the cutoff. The short range repulsiveforce is then F R ij = − ˆ x ij dd ˜ x U sr ij ( x ij ; σ ij , ξ ij ) , (43)where σ ij is the average of the diameters of particles i and j , and in this work we set ξ ij = ξ for all pairs of ions.Note that the repulsion diameter is distinct from the hydrodynamic radius a t discussed in Sec. II E, however, both arerelated to the solvation layer formed around an ion [48]. The potential and repulsive force are illustrated in Fig. 1. IV. DETAILS OF NUMERICAL METHODOLOGY
In this section we discuss the details of the discretization and solution methods of the hydrodynamic and electro-static equations, Eqs. (13) and (32), and the associated particle time stepping algorithm. These are all implementedin the AMReX framework described in Ref. 49 and available at https://amrex-codes.github.io/ . Note that the hydrodynamic equations are in the quasi-hydrostatic approximation (neglect inertial terms) while the electrostatic equa-tions are in the quasi-electrostatic approximation (neglect magnetic fields).
684 0.6 0.7 1.00.90.8 1.1 0.6 0.7 1.00.90.8 1.1432 2100
FIG. 1. Illustration of short range potential and repulsive force due to Eq. (42). For ease of viewing, here we have set ˜ x m = 0 . σ .Note that the actual value used in the simulations, discussed in Sec. V, is ˜ x m = 0 . σ . A. Peskin Kernels
In sections II B and III, we define kernel functions that describe the interaction of the particles with continuousfields; see Eqs. (15), (16), (35), and (36). These functions are represented numerically by Peskin kernels [25], whichapproximate Gaussian functions on a discrete grid. The four-point kernel is defined using δ Pe ( ζ k ) = − | ζ k | + (cid:113) | ζ k | − | ζ k | r , ≤ | ζ k | ≤ − | ζ k | − (cid:113) − | ζ k | − | ζ k | r , < | ζ k | ≤ , < | ζ k | , (44)where ζ = ( x i − r ) / ∆ r , ζ k indicates a single Cartesian component of ζ , ( ζ x , ζ y , or ζ z ), and ∆ r is the grid spacing,which is the same in each direction. Peskin kernels attempt to maximize certain ideal properties of a continuumGaussian, notably isotropy and translational invariance; a detailed discussion is given in Ref. 25. Recently kernelswith a greater number of points in their support, including five and six grid points [50], have been computed. Thesehave improved smoothness and come closer to the ideal case at the cost of increased computation time. Conversely,using a smaller number of points will result in a less accurate kernel but reduced computational cost.Although it is not necessary to use the same kernels for the hydrodynamic and electrostatic fields, for simplicity inthis article we use the same four point kernel in both cases, δ hy ( ζ ) = δ es ( ζ ) = δ Pe ( ζ x ) δ Pe ( ζ y ) δ Pe ( ζ z ) . (45)We further use the same cell size, ∆ r for all fields, although again this is not a necessity.The discretized interpolation operation for a single Cartesian component of the hydrodynamic velocity ( v k ) is J hy h,k ( x i ) v k = ∆ V (cid:88) j ∈ Ω Pe i δ hy ( x i − r j ) v k ( r j ) , (46)where ∆ V = (∆ r ) , and j ∈ Ω Pe i indicates summation over discrete grid points r j that lie within the support of δ hy ( x i − r j ) centered on x i , and on which the discrete velocity is defined. The spreading operation for a single2component of force is given by S hy h ( x i ) F k = N (cid:88) i =1 (cid:88) j ∈ Ω Pe i δ hy ( x i − r j ) F i,k . (47)Similarly, the discretized electrostatic interpolation and spreading operators are, respectively, J es h,k ( x i ) E P k = ∆ V (cid:88) j ∈ Ω Pe i δ es ( x i − r j ) E P k ( r j ) , (48)and S es h ( x i ) q = N (cid:88) i =1 (cid:88) j ∈ Ω Pe i δ es ( x i − r j ) q i . (49)As mentioned above, one reason for using Peskin kernels is to maximize translational invariance. However asmall degree of variance remains – the four point kernel used here has a position dependent hydrodynamic radiusof (1 . ± . r . Were perfect translational invariance achieved, for unbounded/triply periodic domains thestochastic drift term in the particle equation of motion, Eqs. (17) and (19), would be zero everywhere, as the divergenceof the mobility would be zero [40]. However the slight lack of invariance of the kernels renders the stochastic driftterm necessary even in domains without solid boundaries; this is discussed in detail in Ref. 24. B. Fluid solver
We discretize the fluid velocity on a uniform Cartesian grid with mesh spacing ∆ r . We use a staggered grid systemwith normal velocities and force density defined at cell faces, and pressure defined in cell centers, i.e., a standardmarker-and-cell (MAC) discretization [51]. Multiple discretizations are used for differential operators representing thegradient, divergence, and Laplacian. Each of these operators use centered second-order differences. We use ∇ c → fh torepresent the face-centered gradient of a cell-centered field, ∇ ch to represent the cell-centered gradient of a cell-centeredfield, ∇ f → ch · to represent the cell-centered divergence of a face-centered field, ( ∇ h ) f to represent the face-centeredLaplacian of a face-centered field, and ( ∇ h ) c to represent the cell-centered Laplacian of a cell-centered field. Finally, (cid:101) ∇ h · is used to represent the face-centered divergence of a field defined on control volumes corresponding to the shifted(staggered) velocity grid. This is illustrated in two dimensions in Fig. 2. Using these operators, the discrete fluidequation solved at each time step is (see Eqs. (13) and (25)), − η ( ∇ h ) f v + ∇ c → fh p = S hy h F + (cid:114) ηk B T ∆ t ∆ V (cid:101) ∇ h · (cid:98) Z + (cid:20) S hy h (cid:18) x + ∆ R (cid:99) W (cid:19) − S hy h (cid:18) x − ∆ R (cid:99) W (cid:19)(cid:21) k B T ∆ R (cid:99) W , (50a) ∇ f → ch · v = 0 . (50b)Here, (cid:98) Z represents independent normally distributed random numbers defined such that their divergence can becalculated at the locations where v is defined. This is illustrated in two dimensions in Fig. 2, where (cid:98) Z is stored on cellnodes and centers. In three dimensions (cid:98) Z is stored on cell edges and centers. This is described in detail in Ref. 26.The vector (cid:99) W = (cid:110) (cid:99) W , · · · , (cid:99) W i , · · · , (cid:99) W N (cid:111) consists of independent random vectors (cid:99) W i = (cid:16)(cid:99) W i,x , (cid:99) W i,y , (cid:99) W i,z (cid:17) , whereeach (cid:99) W is an independent Gaussian random number of mean zero and variance one.The second line of Eq. (50a) represents the thermal forcing term given by Eq. (24). As discussed in Ref. 24, itis essentially a finite difference representation of the divergence of the spreading operator. It is implemented byspreading a random force, k B T (cid:99) W / ∆ R , at locations randomly offset from the particle positions, x + ∆ R (cid:99) W /
2, andthe negative of that force at positions x − ∆ R (cid:99) W /
2. The distance ∆ R should be smaller than the length scale overwhich S hy varies, but large enough to avoid issues due to numerical roundoff. In this article we have used the value∆ R = 10 − ∆ r .The system of discrete equations formed by Eq. (50) can be solved by a range of approaches, here we have used apreconditioned Generalized Minimal RESidual (GMRES) method [51, 52].3 FIG. 2. Right: A 2D illustration of the MAC [51] discretization used for Eq. (50). Left: 2D illustration of the (cid:101) ∇ h · Z term inEq. (50a), which is the discretization of the term ∇ r · Z from Eq. (11). Random numbers, Z , are generated at the faces ofcontrol volumes around the points where velocities v are defined so that the divergence of the field Z can be calculated. C. Electrostatic solver
The electrostatic potential and charge density, φ and (cid:37) , are defined on a cell centered grid. The Poisson equation,Eq. (32), is discretized by ( ∇ h ) c φ = − (cid:37)/(cid:15). (51)In this case, the resulting system of equations are solved using the geometric multigrid [53] method.The resulting electric field is also stored at cell centers, and is found by E = − ∇ ch φ. (52)The choice of a cell centered differencing in Eq. (52) has some advantages over other options that are worth noting.Importantly, it can be shown that Eq. (52), together with the IB interpolation and spreading, ensures that, at leastfor periodic domains, Newton’s third law is obeyed. That is, the force on ion i due to ion j is equal and opposite tothat on ion j due to ion i , and momentum is conserved. This implies that there is strictly no self-force of one ionon itself, which is an important property desirable in any P3M method. It may appear more natural for staggeredgrids to compute electric field not at cell centers using ∇ ch , but on cell faces using ∇ c → fh , in particular because( ∇ h ) c = ∇ f → ch · ∇ c → fh . However, this choice does not lead to equal and opposite action/reaction forces and cangenerate spurious self-forces.One thing that is lost with the discretisation employed in Eq. (52) is that the force is no longer a negative gradientof an electrostatic potential. That is, the IB-P3M approach followed here cannot be used for applications where bothforces and energies matter and need to be consistent with each other; in BD or DISCOS we only need forces, andsince the system is isothermal there are no issues regarding energy conservation. It is possible to ensure consistencybetween electrostatic forces and energy if one uses the derivative of the Peskin kernel to interpolate − φ at the particlepositions to obtain the electric field/force. However, this requires using a smoothly-differentiable kernel such as the6pt Peskin kernel [50], and also does not lead to momentum conservation.When applying the close range correction to the electrostatic solution (and the steric repulsion), near linear scalingis maintained by using the neighbor list feature included with AMReX. We thank Charles Peskin for sharing with us an analytical proof of this property. Note that this force is, however, not strictly central. That is, the Weak Law of Action and Reaction (Newton’s Third Law) is obeyedbut not the Strong Law of Action and Reaction, which requires that the forces act along the line joining the particles. Note that in the overdamped limit there is no momentum conservation per se, however, the steady Stokes equation would not be solvablein triply-periodic domains if the forces did not sum to zero. D. Temporal Algorithm
The time stepping scheme employed in DISCOS builds on that of the FIB method with the addition of electrostaticforces and dry diffusion. A time step is then defined by the following four steps:1. The charge density (cid:37) is computed on the grid using the spreading operation defined by Eqs. (36) and (49).Equation (51) is then solved using the geometric multigrid method, and Eq. (52) used to obtain the coarseelectric field, E .2. This electric field is interpolated to the particle locations using the operation defined by Eqs. (35) and (48). Thecorresponding force on the particles is found using Eq. (37). For particles at close range, this force is correctedas per Eqs. (38) and (39). Close range WCA interactions are also calculated in this step, and the total force onthe particle calculated as per Eqs. (31)-(43), including the effect of any external fields.3. The force on the particles, F , is spread to the grid storing the force density using the operations defined inEqs. (16) and (47), recovering f s . The random finite difference term is spread as per the second line of Eq. (50a),yielding the thermal forcing f th . The resulting system, Eq. (50), is solved via the GMRES method to computethe velocity at time step n , v n .4. The fluid velocity v n is interpolated to particle locations using the operations defined in Eqs. (15) and (46)to obtain the “wet” component of the particle velocities, Eq. (17). The temporal discretization of the particlediffusion is then given by the midpoint update scheme x n +1 / ,(cid:63)i = x ni + ∆ t J hy h ( x ni ) v n ) , (53) x n +1 i = x ni + ∆ t (cid:34) J hy h ( x n +1 / ,(cid:63)i ) v n + D dry i k B T F ni + (cid:115) D dry i ∆ t W ni . (54)As discussed in Sec. II B, the purpose of the midpoint update is to incorporate the part of the stochastic driftnot accounted for by the random finite difference force density f th .This algorithm is first order in time, as illustrated in Ref. 24, where higher order schemes are also discussed. V. NUMERICAL RESULTS
In this section we test the DISCOS algorithm by comparison with theoretical results for the radial pair correlationfunction and electrical conductivity. Additionally, we analyze the effect of changing the ratio of wet and dry diffusion.In our numerical tests we model a 1:1 strong electrolyte solution, similar to salt water, with species labeled A and B ,with charges q = q A = − q B = 1 . × − C. The solvent is taken to be water at T = 295 K, viscosity η = 0 .
01 g / (cm s),and permittivity (cid:15) = (cid:15) r (cid:15) , with relative permittivity (cid:15) r = 78 . (cid:15) is the vacuum permittivity. The diffusioncoefficients of the ions in water are D tot A = 1 . × − cm / s and D tot B = 1 . × − cm / s, corresponding tohydrodynamic radii of a t = 0 .
185 nm, and a t = 0 .
162 nm, respectively. In all cases the close range potentialparameters are ξ = 10 − ergs, σ = 0 . x m = 0 . A and B , respectively. This justifies the use of the infinite Schmidt numberapproximation discussed above.For each problem the time step, ∆ t , was selected by successive refinement until a negligible change in the resultwas observed. Note that the ion diffusive time scale, a /D , is on the order of 10 picoseconds. The time steps used forthe simulations below were typically constrained to the order of 0.1 ps by the stiffness of the steric and electrostaticinteractions. We are currently examining how the time step may be increased by varying the parameters of thesepotentials.In each case the system size was selected to avoid significant finite size effects. When considering electrostaticinteractions, the relevant length scale is the Debye length [54, 55], λ D = (cid:115) (cid:15)k B T (cid:80) N s j =1 n j q j , (55)where N s it the number of species, and n j is the number density of species j . This is a measure of how far electrostaticeffects persist before they are screened by clouds of opposite charges. For hydrodynamic effects, it should be noted that5the diffusion coefficient of an isolated particle in a triply-periodic domain has well-known strong finite-size correctionsof order a w /L [43], where L is the length of the cubic domain. This is because the mobility of an isolated particle isdetermined by applying a force on it, and a force monopole interacts with its periodic images with 1 /r hydrodynamicinteractions. However each ion diffuses with its Debye cloud, which experiences an equal but opposite force, so thatan ion interacts with its periodic images hydrodynamically as a force dipole [56], which has a finite size correctionof order ( a w /L ) . We are not aware of any systematic analysis of finite-size effects in the literature; here we havesystematically increased the system size until no significant change in the result is observed. A. Radial Pair Correlation Function
We first validate the equilibrium properties of the DISCOS method by measuring radial pair correlation functionsbetween ions of like and opposite charge. In general, the radial pair correlation function is the normalized time-average density of particles as a function of radius from an arbitrary reference particle. For a binary system, the paircorrelation function between species α and β is defined as g αβ (˜ x ) = lim τ →∞ V N αβ π ˜ x τ (cid:90) τ N (cid:88) i,j,i (cid:54) = j δ (˜ x i − x ij ) δ α,s i δ β,s j dt, (56)where V is the system volume, s i and s j are the species of particles i and j , N αβ is the number of species pairs, and˜ x i is the radial distance from particle i . In practice, for each snapshot in time, for each ion of type α , we count thenumber of ions of type β in thin spherical shells with a specified bin width, and compute a number density in eachbin by dividing by the volume of the shell. We average the results over all α ions and over many time steps, andnormalize the result with the average number density of β .For low to moderate ion concentrations Debye-H¨uckel theory [54, 55] gives g αβ (˜ x ) ≈ exp( − U αβ (˜ x ) /k B T ) , (57)where the potential is U αβ (˜ x ) = q α q β π(cid:15) e − ˜ x/λ D ˜ x + U sr αβ (˜ x ) . (58)The first term in Eq. (58) is the screened Coulomb potential [54, 55] and the second is the short-ranged repulsion.In Figs. 3 and 4, we compare the approximate theoretical expression of Eq. (57) with results obtained from DISCOS.In Fig. 3, the comparison is shown for molarity (moles of cation or anion per liter of solvent) of 0.1 M; excellentagreement is observed. Fig. 4 shows the comparison for a molarity of 1.0 M for two different ratios of wet and drydiffusion. A negligible difference is observed between second and third simulations; this is unsurprising as we cansee from Eq. (57) that the pair correlation function does not depend on the hydrodynamic properties of the solvent.Reasonable agreement is observed, with the largest deviation observed in the peak of the opposite charge result.We note that Eq. (57) is derived in the low concentration limit so we expect decreasing agreement with increasingconcentration.At the higher concentrations the Percus-Yevick and hypernetted chain (HNC) approximations are moreaccurate [57]. B. Electrical Conductivity
For an applied electric field of magnitude E the total current density in an electrolyte solution is I = N s (cid:88) j =1 n j q j µ j E, (59)where µ j is the mobility of species j . The conductivity is defined by Ohm’s law as C = I /E and from the Nernst-Einstein relation [55] the mobility in the limit of infinite dilution is µ j = D tot j /k B T . In this limit, the uncorrectedtotal conductivity C is the sum of the species conductivities C i , C = N s (cid:88) j =1 C j = N s (cid:88) j =1 n j q j D tot j k B T . (60)6
FIG. 3. Radial distribution for a molarity of 0.1 M. The right plot shows the pair correlation for ions of like charge; alsoindicated is the Debye length. The left plot shows the pair correlation for ions of opposite charge. The solid line shows theapproximate analytical solution given by Eq. (57), and the squares show the numerical result from DISCOS. These results werecollected using a bin width of 0.025 nm, with a sample size of 50,000 – the simulation was run for 100,000 steps, with the first50,000 steps used for equilibration. For visual clarity, only every fourth point is displayed. The grid size was set such that thetotal diffusion of species A was 14.7% wet and 85.3% dry, and for species B , 12.9% wet and 87.1% dry. This corresponds to∆ r = 1nm, with a 32 × ×
32 cell periodic domain, and 3946 ions. The time step was ∆ t = 0 . Debye-H¨uckel-Onsager theory gives the electrolyte conductivity [55, 58] as the sum of the uncorrected total conduc-tivity, electrophoretic contribution, and a contribution due to the relaxation effect, C = C + C ep + C relx . (61)See Ref. 12 for the relationship of this theory to fluctuating hydrodynamics.The electrophoretic contribution is due to screening charges of opposite sign about each ion imparting a retardingviscous stress on the ion, given by C ep = − N s (cid:88) j =1 n j q j πηρλ D λ D λ D + a DH j . (62)The last term is a Debye-H¨uckel correction to account for the finite ion size, which we take to be a DH j = 0 . C relx = − (cid:88) j =1 C π (2 + √ q j (cid:15)k B T λ D λ D λ D + a DH j . (63)For a DH j (cid:28) λ D both the electrophoretic and relaxation contributions go as the square root of the ionic concentration.The C term contains no hydrodynamic interactions between ions; it is therefore correctly captured by a simulationwith any wet/dry diffusion ratio. The same is true for the C relx term, as this arises entirely from electrostaticinteractions which are well-resolved using the P3M approach. The C ep term arises from the collective hydrodynamiceffect of an ion and its screening cloud of opposite charge, it is therefore only completely captured by a simulationusing 100% wet diffusion; the relative importance of this effect is discussed in the next section.7 FIG. 4. Radial distribution for a molarity of 1.0M. The right plot shows the pair correlation for ions of like charge; alsoindicated is the Debye length. The left plot shows the pair correlation for ions of opposite charge. The solid line shows theapproximate analytical solution given by Eq. (57), and the boxes and circles show numerical result from DISCOS. The numericalparameters are the same as those described for Fig. 3, except that two different grid resolutions are shown. The black circlesshow the case where the grid size was set such that the total diffusion of species A was 14.7% wet and 85.3% dry, and forspecies B , 12.9% wet and 87.1% dry. This corresponds to ∆ r = 1nm, with a 16 × ×
16 cell periodic domain, and 4932 ions.The red squares show the case where the grid size was set such that the total diffusion of species A was 58.9% wet and 41.1%dry, and for species B , 51.8% wet and 49.2% dry. This corresponds to ∆ r = 0 . × ×
32 cell periodic domain,and 618 ions.
When the applied electric field is sufficiently strong the distortion of the ionic clouds results in an increased electricalconductivity, which is known as the first Wien effect. For a 1:1 strong electrolyte solution the Wien corrections tothe conductivity are [59] C relxW ( a ) = A W ( a ) C relx , (64) C epW ( a ) = B W ( a ) C ep , (65)where A W ( a ) = 3(1 + (cid:112) / a (cid:16) a (cid:112) a − √ a + arctan( √ a ) − arctan( a / (cid:112) a ) (cid:17) , (66) B W ( a ) = 1 √ a (cid:16) a arcsinh( a ) − a (cid:112) a + √ a − (1 + 2 a ) (cid:104) arctan( √ a ) − arctan( a / (cid:112) a ) (cid:105) (cid:17) , (67)and a = ( qλ D E ) / ( k B T ) is a dimensionless parameter for the electric field strength. These correction factors are shownin Fig. 5. Note that A W ( a ) → B W ( a ) → / √ a → ∞ .In the presence of an applied electric field, the conductivity vector can be computed by C = lim τ →∞ E V τ [ Z ( τ ) − Z (0)] , (68)where the ion polarization (displacement of the “center of charge”) is Z ( t ) = N (cid:88) i =1 q i x i ( t ) . (69) Note that there is a “second Wien effect” which applies to weak electrolytes. FIG. 5. Wien corrections to the DHO relaxation effect and electrophoretic contribution as a function of a = ( qλ D E ) / ( k B T ). An alternative way to measure conductivity in the absence of an applied field is given by the Einstein-Helfandformula [60], in terms of the long-time diffusion coefficient of the center of charge, C EH = lim τ →∞ k B T V τ (cid:90) τ [ Z ( t ) − Z (0)] dt. (70)This is derived from the Green-Kubo relation [61, 62], and is discussed in Refs. 60 and 63.Formally, Eqs. (68) and (70) hold when τ is large enough to ensure that the short time correlations have decayed;increasing τ beyond this point has no effect on the statistical convergence, as the variance of Z does not decrease withtime. An accurate result can be obtained by increasing N , or averaging an ensemble of simulations. We find that themean-square displacement of the center of charge is perfectly linear in time on our overdamped time scales for all τ examined, so the conductivity is not sensitive to the choice of τ . The correct result can therefore be obtained by setting τ ∼ ∆ t to maximize the statistical sample size. To ensure de-correlation between measurements the conductivity wassampled at intervals of 10 ps, which is on the order of the ion diffusion time, when calculating the time average inEq. (68). For Eq. (70), the displacement in Z was measured and then reset every 10 ps. Note that this has no impacton the result itself, but it influences the calculation of the statistical error.Three sets of conductivity measurements were performed at a concentration of 0.1 M: (i) Measurements with zeroapplied field using Eq. (70), (ii) using a weak applied field of 10 V / cm ( a = 0 . V / cm ( a = 37 . D wet and D dry , were tested. Tests using case (i) with Eq. (70) were also performed at concentrations of 0.0115 M,0.5 M, and 2 M. The results are summarized in Table I and illustrated in Fig. 6.All simulations were performed in a cubic domain with length L = 10 .
043 nm. For concentrations 0.0115 M, 0.1 M,0.5 M, and 2 M, the domain contained 14, 122, 610, and 2440 ions, respectively. For comparison, a molecular dynamicssimulation of this system would contain roughly 30,000 water molecules.To vary the wetness, the computational domain was divided into 8 , , , or 64 grid cells for both the Stokesand Poisson solvers; simulations with fully dry diffusion used a 32 grid for the Poisson solver. The 8 grid cell casescorrespond to wetness percentages of 11 .
8% and 10 .
3% for ion A and B, respectively, while the 64 grid cell casescorrespond to wetness percentages of 93 .
8% and 82 . t was typically between 0.01 and 0.02 ps.It is important to note that while a Stokes grid as coarse as 8 does little to accelerate the far-field calculations fornearest image particles, the use of a periodic Stokes solver is crucial to account for the hydrodynamic interactionswith periodic image particles; the same applies also to the Poisson solver.In all cases, good agreement is observed between the results with the highest D wet value and the theoreticalprediction, for both high and low field cases, and at all concentrations. The weak and strong field results are inparticularly good agreement with the theory for both fully dry and fully wet simulations. Note that DHO theory may9not be accurate at high molarities, especially since we do not tune the value of the finite-size radius a DH for eitherspecies. The dry ( D wet =0) cases also compare reasonably well with the theoretical result with the electrophoreticcomponent removed, in agreement with the discussion above.For the weak field simulations, the strength of the field has been set such that the Wien effect is negligible; wetherefore expect the same results for the weak and zero field cases. Interestingly, a difference of about 3% is observedfor the highest wet percentage, and 5% for the lowest. The cause of this discrepancy is not apparent, but warrantsfurther investigation. Measured. Wet% A , Wet% B Theory 0% 11.75% 23.5% 46.9% 93.8% TheoryNo EP 0% 10.33% 20.65% 41.3% 82.6% With EPZero E -field0.01147 M 0.106 0.1086(18) 0.1013(6) 0.1020(21) 0.1002(21) 0.1010(18) 0.1000.1 M 0.897 0.946(12) 0.869(27) 0.837(21) 0.813(9) 0.803(12) 0.7760.5 M 4.33 4.70(6) 4.21(9) 3.95(6) 3.63(6) 3.54(15) 3.342 M 16.8 18.5(6) 16.7(9) 14.8(9) 12.8(6) 11.6(9) 11.4Weak E -field0.1 M 0.897 0.898(6) 0.831(6) 0.785(6) 0.785(9) 0.777(6) 0.776Strong E -field0.1 M 0.943 0.943(3) 0.870(3) 0.851(3) 0.841(3) 0.837(3) 0.850TABLE I. Conductivity measurements in siemens per meter for zero E -field (by Einstein-Helfand, Eq. (70)) and weak and strong E -field (center of charge displacement, Eq. (68)). The uncertainty notation is defined such that the numbers in parenthesisindicates the error bar for the final digits, e.g., 0 . . ± . . . ± . E -field cases. For comparison with the0% wet case, the theory prediction without the electrophoretic (EP) correction is also included in the table. A plot of theseresults is given in Fig. 6. C. Effect of wet vs. dry diffusion
The results in the previous section show that while some wet diffusion is needed to capture the electrophoreticcontribution to conductivity, in most cases increasing the wet diffusion above 50% has a relatively small impact. Fur-thermore the relative effect depends on concentration, with wet diffusion making a greater contribution at higher molarconcentrations; this can be clearly seen in Fig. 6. As mentioned in Sec. II E, reducing the hydrodynamic resolution(by reducing D wet ) neglects short range hydrodynamic interactions. It seems clear that short range hydrodynamicinteractions should increase in importance as concentration increases, i.e., as the particles are on average closer to-gether. However it must be noted that this will also increase the relative importance of electrostatic interactions,which decay more rapidly than hydrodynamic influence.To estimate the importance of these effects, and obtain an estimate for the required D wet value for a given con-centration, we consider two ions, 1 and 2. The electrostatic force between them is given by the Coulomb interaction,Eq. (38), while the hydrodynamic interaction is given by the RPY tensor, Eq. (7). To represent the dry componentof the particle mobility, we introduce the dry tensor M dry ( a ) = 16 πηa d I . (71)The velocity of particle 1 has two contributions: That directly induced by the Coulomb force, V dir1 = (cid:16) M RPY (0 , a w ) + M dry ( a d ) (cid:17) F C12 = M RPY (0 , a t ) F C12 = M dry ( a t ) F C12 . (72)and the disturbance in the fluid induced by particle 2, which is subject to the opposite force, V dst1 = − M RPY ( x , a w ) F C12 , (73)with the total velocity of particle 1 (parallel to the vector x ) being V = V dir1 + V dst1 . (74)0 FIG. 6. Relative difference in conductivity from DHO theory, as a function of wetness percentage, for a range of molarities;based on data from Table I.
In Fig. 7 we plot V as a function of particle separation. Five curves are shown, ranging from 0% to 100% wet,with increments of 25%. Also marked, with vertical lines, are the measured average nearest particle separations fortwo concentrations, 0.1 M (black) and 2 M (red). At 0.1 M varying the wetness from 50-100% produces only a 2%change in particle velocity. At 2 M the same range produces a change in velocity of 33%, suggesting a higher wetnesspercentage is necessary to accurately capture the solution at higher molarities, as expected. FIG. 7. Particle velocity due to electrostatic and hydrodynamic interactions between two ions, as a function of separation.The velocity V is calculated using Eqs. (72)-(74), using a radius averaged from species A and B , i.e., ¯ a t = ( a At + a Bt ) / a w = ( a Aw + a Bw ) /
2, ¯ a d = ( a Ad + a Bd ) /
2. A range of wet percentages are shown, the corresponding particle radii can be found withEqs. (29) and (30) and the radii quoted at the start of Section V. The measured average closest ion separations correspondingto concentrations 0.1 M and 2 M are indicated with the black and red vertical lines. x = Γ( )( πn ) − / ≈ . n − / , (75)where Γ( x ) is the gamma function. Assuming a cubic lattice arrangement gives x = n − / . (76)Also shown are the conductivity measurements from Table I, again normalized by the 0% wet value. For both 0.1 Mand 2 M the conductivity measurements follow a similar trend to the particle velocity , suggesting that the simpleparticle velocity calculation from Eqs. (72)-(74) is a useful indicator of the lower bound of wet diffusion necessary toaccurately simulate a given molar concentration. For example, the change of 2% in the particle velocity for the 0.1 Mcase corresponds to a 1% change in conductivity for the two highest wetness ratios, while the 33% change in velocityfor the 2 M case corresponds to a change in conductivity of 10%. FIG. 8. Normalized particle velocity due to Eq. (74) (lines), and normalized conductivity measurements from Table I (symbols),as a function of wet/dry ratio of species A. All values are normalized by the 0% wet value. Black signifies 0.1 M, red 2 M. Thecurves indicate the particle velocity, calculated at the average nearest particle separation measured from simulation (solid), andcorresponding to random (dashed), and cubic lattice (dotted) configurations. Squares are zero-field conductivity measurementsusing the Einstein-Helfand formula (Eq. (70)), diamonds indicate weak-field conductivity measurements using the center ofcharge displacement (Eq. (68)), and crosses are the prediction of DHO theory (Eq. (61)).
VI. CONCLUSIONS
In this article we have presented a new approach for the simulation of electrolytes at the mesoscale. The methodextends the FIB method to include electrostatic interactions, and incorporates a “dry diffusion” process to correctfor under-resolution of the hydrodynamics. The dry diffusion term allows the method to be applied to systems ofelectrolytes where application of the standard FIB approach would require an excessively fine grid to resolve particlesof nanoscopic size. We note that our approach is similar to the Stochastic Eulerian Lagrangian Method [27, 65, 66]and the stochastic Force Coupling Method [67–69]; the relationship of each of these to the FIB method is discussedin Ref. 24. In principle, some of the techniques outlined in this article could also be applied in the context of thesedistinct numerical approaches.2In sections V A and V B, we have demonstrated that the DISCOS method accurately reproduces several importantproperties of electrolytes. In Sec. V C, we have demonstrated that the dry diffusion approach can effectively replaceresolving the near-field hydrodynamics at scales smaller than the typical ion-ion distance. At 0.1 M, reducing thewetness to approximately 50% produces only a 2% change in the conductivity, but at 2 M a 10% change is observed.Our results clearly demonstrate that, for denser electrolyte solutions, near-field hydrodynamic interactions con-tribute to macroscopic transport properties such as conductivity. As discussed in great detail in Ref. 40, wet diffusion,which is diffusion due to advection by a fluctuating velocity field, and standard dry diffusion, are very different bothphysically and mathematically. This difference cannot be seen for a single particle since ultimately only the total(effective or renormalized) diffusion coefficient matters; it is worth recalling that in Ref. 29 only a single nano-particleis analyzed. The difference between wet and dry diffusion, however, manifests itself in collective properties such asthe spectrum of nonequilibrium concentrations [40]. However, experimentally distinguishing between wet and drydiffusion based on this spectrum requires examining length scales below light-scattering reach, and, at the scales ob-servable in typical experiments, one only sees a total (effective or renormalized) diffusion coefficient. Our work showsthat electrical conductivity, which is easily measured to high accuracy, is a sensitive probe that is affected by thehydrodynamic interactions (equivalently, the spectrum of the fluctuating velocity field [40]) at distances comparableto the typical ion-ion spacing. In the future, it is important to perform more detailed comparisons between simula-tions for different degrees of “wetness” and results for denser electrolytes obtained using experiments or large-scalemolecular dynamics simulations. This will shed light on how good of a model the Rotne-Prager-Yamakawa tensor isfor solvated ions.A key advantage of the FIB method is the ability to handle other types of boundary conditions, notably confinementby no-slip walls. As the discussion of DISCOS in this context requires considerable exposition, of both details of thealgorithm and the physics that arises, we defer this to a future publication. It is perhaps more relevant to contrastthe numerical method used in DISCOS to the state-of-the-art Positively Split Ewald (PSE) method for BrownianDynamics with hydrodynamic interactions developed in Ref. 21 by using some ideas from spectral Ewald methodsfor electrostatics [46] combined with the idea of using fluctuating hydrodynamics to generate the far-field Browniandisplacements. The PSE method includes a spectrally-accurate implementation of the Hasimoto-Ewald splitting usedin GGEM [17]; this splitting allows one to choose the grid size arbitrarily by using the same idea of local near-fieldcorrections that we used here for the Poisson equation. However, at present, the PSE approach relies on Fourier-spacedecompositions and using the FFT algorithm, and therefore only apply to triply-periodic domains. It remains afuture challenge to find a way to incorporate near-field hydrodynamic corrections to correct for under-resolution ofthe hydrodynamics by the Stokes solver used in FIB and DISCOS, while also computing Brownian motion withoutexpensive iterations.In future work we intend to apply DISCOS to strong electrolyte solutions that are acids (e.g., HCl) or bases (e.g.,NaOH) instead of salts. Also of interest are solutions in which the Wien effect is more pronounced, such as weakelectrolyte acids and bases (e.g., acetic acid, ammonium hydroxide) and strong electrolyte salts that are not 1:1 (e.g.,MgCl ). Furthermore, recent work suggests that non-aqueous solvents with low permittivity (e.g., benzene) exhibitunusual Wien effects [70].Here we have described the DISCOS approach for simple mono-atomic, 1:1 electrolytes, where each ion is describedby a single IB kernel. In principle, multiple kernels may be combined to form complex molecules such as chargedpolymers or DNA strands. This approach has already been tested in the context of BD, for example in Refs. 30 and31, and for electrically neutral systems using FIB, for example in Refs. 71–73. Future work may also include theextension of DISCOS to encompass these cases. Further, there is the possibility of simulating large molecules, suchas DNA, in a solution represented as a continuum that includes smaller ionic species. This was done previously usinga deterministic solution in Ref. 74, in this case the solution would be modeled using the fluctuating hydrodynamicapproaches outlined in Refs. 11–13. ACKNOWLEDGEMENTS
The authors thank Changho Kim for discussions regarding electrolyte chemistry, and Charles Peskin and JonathanGoodman for discussions regarding the IB-P3M method. This work was supported by the U.S. Department of Energy,Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under contractNo. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center,a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy underContract No. DE-AC02-05CH11231. A. Donev was supported in part by the Division of Chemical, Bioengineering,Environmental and Transport Systems of the National Science Foundation under award CBET-1804940.3
Appendix A: Accuracy of immersed-boundary P3M solver
In this appendix we examine the accuracy of the immersed-boundary P3M method described in Sections III Aand IV C. In particular, we focus on the degree of translational and rotational invariance of the electrostatic forcescomputed by the method for a pair of ions of equal and opposite charge q i = − q j at a distance r from each other,using a uniform regular grid of spacing ∆ r . The periodic box size L (cid:29) r so that periodic effects can be neglected.We include in our tests both the 4-point (see Eq. (44)) and 6-point [50] Peskin kernels. In order to compare toalternative methods, these tests were performed in Matlab using FFTs instead of multigrid to solve the (discrete)Poisson equation. In addition to the 2nd-order electrostatic solver based on the 7-point Laplacian ( ∇ h ) c and centereddifference ∇ ch , we also include, for comparison, results obtained by using a fourth-order compact isotropic 24-pointLaplacian in Eq. (51) and a 4-point centered difference in Eq. (52).In addition to the fourth-order finite-difference Poisson solver, we also implement a spectral Fourier solver, as usedin the Spectral Ewald (SE) method [46]. As explained in Section III A, the SE approach (and essentially all otherEwald methods) uses a Gaussian of standard deviation σ as the kernel δ es in Eqs. (35) and (36). This allows oneto compute the electrostatic force between the two Gaussian charge clouds analytically for an unbounded domain, giving the required local correction F LCGaussian = (cid:18) ∆ rr (cid:19) (cid:18) erfc (cid:16) r σ (cid:17) + rσ √ π exp (cid:18) − r σ (cid:19)(cid:19) (A1)that decays exponentially with distance. However, achieving 3-4 digits of accuracy with Gaussian kernels requiresthe support of the kernel to be on the order of 10-12 grid points in each dimension, which is expensive. Therefore,to make a more fair comparison with the 6-point Peskin kernel, here we couple the spectral Poisson solver with arecently-developed “exponential of a semicircle” (ES) kernel [75] with width w = 6; we optimized the value of theparameter β = 1 . w for translational and rotational invariance. Figure 9 shows numerical data for the radial component of the near-field correction to the electrostatic force F LC ij for several kernels and Poisson solvers, expressed in terms of the dimensionless function F LC ( x = r/ ∆ r ) appearing inEq. (41). To quantify the loss of translational and rotational invariance, we randomly generate the position of manypairs of particles at a distance r from each other. Note that if F LC ij were purely radial and its magnitude were purelya function of r , the local-correction approach used here could completely correct the electrostatic force computed bythe discrete solver to match the Coulomb force between point charges. Therefore, the main cause of inaccuracy isthe loss of translational and rotational invariance. In Fig. 9 we also show the correction F LC ( x = r/ ∆ r ) for a pairof Gaussian charge clouds given in Eq. (A1). The value of the Gaussian width σ was fitted to match the numericalresults over the range 0 < x <
3, giving σ ≈ . r for the 4-point and σ ≈ . r for the 6-point Peskin kernelswith the 2nd order Poisson solver. This can be thought of as an effective width of the charges similar to the effectivehydrodynamic radius a w for Stokes flow.The spread in the symbols in Fig. 9 shows the loss of translational and rotational invariance. We see in the toppanel that the mean correction force decays algebraically rather than exponentially as it does for Gaussian charges.We also see that the 6pt kernel improves the invariance, especially at short distances. Using a fourth-order isotropicPoisson solver further improves the invariance for all distances (data not shown in Fig. 9), as does using a spectralPoisson solver with the 6-point ES kernel, as shown in the bottom panel of Fig. 9. For the 4-point kernel, we see inthe top panel of the figure that truncating the local corrections for pairs of particles further than 3-grid points apart, ψ = 3∆ r , makes the error due to loss of invariance dominate the error due to truncation of F LC ( x = r/ ∆ r ), andtherefore there is no point in further increasing the cutoff distance ψ ; similarly, for the 6-point kernel ψ = 4∆ r issuitable.We further quantify the amount of translational and rotational variance in Fig. 10. The top panel shows the spreadin the symbols in Fig. 9 for each distance r , measured as two standard deviations of the samples. We normalize thespread by 1 /x to obtain a measure of the invariance of the radial component of the electrostatic force relative tothe true Coulomb potential, and express the error in percent; note that this is the error in the final answer for theelectrostatic force between two ions. The bottom panel shows the magnitude of the non-radial component of theelectrostatic force computed by the Poisson solver, F P ij , relative to it’s radial component, again expressed as percenterror. The top panel of Fig. 10 shows that the second-order Poisson solver is translationally and rotationally invariantto about 8% for the 4-point kernel and about 6% for the 6-point Peskin kernel, with the 6-point kernel being notably Note that the periodic boundaries are implemented by using a Fourier series instead of a Fourier sum. We point out again that this approach, while similar to the Spectral Ewald method [46] based on the non-uniform FFT [75] using theES kernel, is different because we do use the same kernel to both smear the charges and to communicate to the grid, rather than usingtwo different kernel functions (Gaussian for smearing and ES for grid transfers in the NUFFT). x = r/ ∆ r -5 -4 -3 -2 -1 C o rr ec ti on F L C ( x ) σ =0.78 ∆ r x = r/ ∆ r -5 -4 -3 -2 -1 C o rr ec ti on F L C ( x ) σ =0.89 ∆ r FIG. 9. Local correction function F LC ( x = r/ ∆ r ) computed from the radial component of the near-field correction to theelectrostatic force, F LC ij , for the Peskin 4pt kernel (top panel) and the Peskin 6pt kernel (bottom panel). Red circles aresamples for different positions and orientations of the pair of ions, while dashed blue lines indicate the continuum theoreticalnear-field correction F LCGaussian for Gaussian charges of standard deviation σ indicated in the legend (see Eq. (A1)). The blacksolid line in the top panel indicates the empirical mean used to tabulate F LC in our method. The bottom panel includesresults (black squares) for a spectral Poisson solver with a 6-point exponential of a semicircle (ES) kernel [75] with parameter β = 1 . × more accurate at short distances. The fact that errors become comparable for large distances regardless of the widthof the kernel suggests that for large distances the error is dominated by the anisotropy of the discrete Green’s functionfor the standard 7-point Laplacian. Switching to a fourth-order Laplacian drops the maximum error down to 2% forthe 6-point kernel, while the spectral solver with the 6-point ES kernel drops the error below 0.5%. The bottom panelof the figure shows that all combinations of kernels and solvers studied here achieve a central electrostatic force towithin 1.5%, which is remarkably good. r / ∆ r % e rr o r du e t o i nv a r i a n ce r / ∆ r % ce n t r a lit y FIG. 10. Percent error due to loss of translational and rotational invariance in the immersed-boundary P3M method based ondata shown in Fig. 9. In addition to the method used in this work, we include for comparison the improvement gained by usinga fourth-order isotropic Poisson solver, as well as a spectral Poisson solver with a 6-point ES kernel. (Top) Percent error in theradial component of the total electrostatic force due to loss of invariance. (Bottom) Magnitude of the non-radial component ofthe IB force F P expressed as percent of the radial component. The symbols show the empirical mean over samples of pairs ofpoints, while the dashed lines show the maximum value over the samples. [1] A. J. Grodzinsky, Fields, forces, and flows in biological systems (Garland Science, 2011).[2] B. J. Kirby,
Micro and nanoscale fluid mechanics: Transport in microfluidic devices (Cambridge university press, 2010).[3] E. Verpoorte, Electrophoresis , 677 (2002).[4] L. Li, X. Wang, Q. Pu, and S. Liu, Analytica Chimica Acta , 1 (2019).[5] J. C. Bachman, S. Muy, A. Grimaud, H.-H. Chang, N. Pour, S. F. Lux, O. Paschos, F. Maglia, S. Lupart, P. Lamp, et al. ,Chemical reviews , 140 (2016).[6] K. Scott, Renewable and Sustainable Energy Reviews , 1406 (2018).[7] M. Andersson, J. Yuan, and B. Sund´en, Applied Energy , 1461 (2010).[8] T. Jahnke, G. Futter, A. Latz, T. Malkow, G. Papakonstantinou, G. Tsotridis, P. Schott, M. Grard, M. Quinaud,M. Quiroga, A. Franco, K. Malek, F. Calle-Vallejo, R. F. de Morais], T. Kerber, P. Sautet, D. Loffreda, S. Strahl,M. Serra, P. Polverino, C. Pianese, M. Mayur, W. Bessler, and C. Kompis, Journal of Power Sources , 207 (2016).[9] A. Arsalis, Renewable and Sustainable Energy Reviews , 391 (2019).[10] L. Landau and E. Lifshitz, Course of theoretical physics. Vol. 6: Fluid mechanics (London, 1959).[11] J.-P. P´eraud, A. Nonaka, A. Chaudhri, J. B. Bell, A. Donev, and A. L. Garcia, Physical Review Fluids , 074103 (2016).[12] A. Donev, A. L. Garcia, J.-P. P´eraud, A. J. Nonaka, and J. B. Bell, Current Opinion in Electrochemistry , 1 (2019).[13] A. Donev, A. J. Nonaka, C. Kim, A. L. Garcia, and J. B. Bell, Physical Review Fluids , 043701 (2019).[14] D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications , Vol. 1 (Elsevier, 2001).[15] J. Bockris, M. Devanathan, and K. M¨uller, in
Electrochemistry (Elsevier, 1965) pp. 832–863.[16] M. Fixman, The Journal of Chemical Physics , 1527 (1978).[17] J. P. Hern´andez-Ortiz, J. J. de Pablo, and M. D. Graham, Physical review letters , 140602 (2007).[18] X. Zhao, J. Li, X. Jiang, D. Karpeev, O. Heinonen, B. Smith, J. P. Hernandez-Ortiz, and J. J. de Pablo, The Journal ofChemical Physics , 244114 (2017).[19] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (Taylor & Francis, Inc., Bristol, PA, USA,1988).[20] S. Srinivasan and A.-K. Tornberg, Research in the Mathematical Sciences , 35 (2018).[21] A. M. Fiore, F. B. Usabiaga, A. Donev, and J. W. Swan, J. Chem. Phys. , 124116 (2017), software available at https://github.com/stochasticHydroTools/PSE .[22] T. Ando, E. Chow, Y. Saad, and J. Skolnick, The Journal of Chemical Physics , 064106 (2012).[23] E. Chow and Y. Saad, SIAM Journal on Scientific Computing , A588 (2014).[24] S. Delong, F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev, The Journal of Chemical Physics ,134110 (2014).[25] C. S. Peskin, Acta Numerica , 479 (2002).[26] A. Nonaka, Y. Sun, J. Bell, and A. Donev, Communications in Applied Mathematics and Computational Science , 163(2015).[27] P. J. Atzberger, Journal of Computational Physics , 2821 (2011).[28] L. af Klinteberg, D. S. Shamshirgar, and A.-K. Tornberg, Research in the Mathematical Sciences , 1 (2017).[29] P. Espa˜nol and A. Donev, The Journal of Chemical Physics , 234104 (2015).[30] R. M. Jendrejack, J. J. de Pablo, and M. D. Graham, The Journal of Chemical Physics , 7752 (2002).[31] R. M. Jendrejack, D. C. Schwartz, M. D. Graham, and J. J. de Pablo, The Journal of Chemical Physics , 1165 (2003).[32] F. Balboa Usabiaga, X. Xie, R. Delgado-Buscalioni, and A. Donev, The Journal of Chemical Physics , 214113 (2013).[33] A. M. Fiore and J. W. Swan, Journal of Fluid Mechanics , 544 (2019).[34] Y. L. Klimontovich, Physica A: Statistical Mechanics and its Applications , 515 (1990).[35] M. H¨utter and H. C. ¨Ottinger, Journal of the Chemical Society, Faraday Transactions , 1403 (1998).[36] R. Kekre, J. E. Butler, and A. J. Ladd, Physical Review E , 011802 (2010).[37] J. Rotne and S. Prager, The Journal of Chemical Physics , 4831 (1969).[38] M. Fixman, Macromolecules , 1204 (1986).[39] P. S. Grassia, E. J. Hinch, and L. C. Nitsche, Journal of Fluid Mechanics , 373 (1995).[40] A. Donev, T. G. Fai, and E. Vanden-Eijnden, Journal of Statistical Mechanics: Theory and Experiment , P04004(2014).[41] G. Da Prato, Kolmogorov Equations for Stochastic PDEs (Birkhauser, 2004).[42] A. Donev, E. Vanden-Eijnden, A. L. Garcia, and J. B. Bell, Comm. Appl. Math and Comp. Sci. , 149 (2010).[43] F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev, Comput. Methods Appl. Mech. Engrg. , 139(2014), code available at https://github.com/fbusabiaga/fluam .[44] P. Ahlrichs and B. D¨unweg, The Journal of Chemical Physics , 8225 (1999).[45] P. Ahlrichs and B. D¨unweg, International Journal of Modern Physics C , 1429 (1998).[46] D. Lindbo and A.-K. Tornberg, Journal of Computational Physics , 8744 (2011).[47] J. D. Weeks, D. Chandler, and H. C. Andersen, The Journal of Chemical Physics , 5237 (1971).[48] R. H. Stokes and R. A. Robinson, Journal of the American Chemical Society , 1870 (1948).[49] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, M. Katz,A. Myers, T. Nguyen, A. Nonaka, M. Rosso, S. Williams, and M. Zingale, Journal of Open Source Software , 1370 (2019). [50] Y. Bao, J. Kaye, and C. S. Peskin, Journal of Computational Physics , 139 (2016), software and updated documentationavailable at https://github.com/stochasticHydroTools/IBMethod , including also a new 5-pt kernel with three continuousderivatives.[51] M. Cai, A. Nonaka, J. B. Bell, B. E. Griffith, and A. Donev, Communications in Computational Physics , 1263 (2014).[52] Y. Saad and M. H. Schultz, SIAM Journal on scientific and statistical computing , 856 (1986).[53] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes 3rd edition: The art of scientificcomputing (Cambridge university press, 2007).[54] P. Debye and E. H¨uckel, Physikalische Zeitschrift , 185 (1923).[55] R. A. Robinson and R. H. Stokes, Electrolyte Solutions: Second Revised Edition , Dover Books on Chemistry Series (DoverPublications, Incorporated, 2012).[56] D. Long and A. Ajdari, The European Physical Journal E , 29 (2001).[57] J. Hansen and I. McDonald, Theory of Simple Liquids (Elsevier Science, 1990).[58] L. Onsager and R. Fuoss, J. Phys. Chem. , 2689 (1932).[59] H. Harned and B. Owen,
The physical chemistry of electrolytic solutions: by Herbert S. Harned and Benton B. Owen ,American Chemical Society. Monograph series, no. 137 (Reinhold Publishing Corporation, 1958).[60] M. Jardat, O. Bernard, P. Turq, and G. R. Kneller, The Journal of Chemical Physics , 7993 (1999).[61] R. Kubo, Reports on progress in physics , 255 (1966).[62] R. E. Jones and K. K. Mandadapu, The Journal of chemical physics , 154102 (2012).[63] K. Nieszporek, J. Nieszporek, and M. Trojak, Computational and Theoretical Chemistry , 52 (2016).[64] W. Feller, An introduction to probability theory and its applications , Vol. 2 (John Wiley & Sons, 2008).[65] P. J. Atzberger, P. R. Kramer, and C. S. Peskin, Journal of Computational Physics , 1255 (2007).[66] P. Plunkett, J. Hu, C. Siefert, and P. J. Atzberger, Journal of Computational Physics , 121 (2014).[67] M. Maxey and B. Patel, International journal of multiphase flow , 1603 (2001).[68] S. Lomholt and M. R. Maxey, Journal of Computational Physics , 381 (2003).[69] E. E. Keaveny, Journal of Computational Physics , 61 (2014).[70] D. Lesnicki, C. Y. Gao, B. Rotenberg, and D. T. Limmer, Phys. Rev. Lett. , 206001 (2020).[71] B. Sprinkle, A. Donev, A. P. S. Bhalla, and N. Patankar, arXiv preprint arXiv:1901.06427 (2019).[72] B. Sprinkle, F. Balboa Usabiaga, N. A. Patankar, and A. Donev, The Journal of Chemical Physics , 244103 (2017).[73] F. Balboa Usabiaga, B. Kallemov, B. Delmotte, A. Bhalla, B. Griffith, and A. Donev, Communications in AppliedMathematics and Computational Science , 217 (2017).[74] J. P. Hern´andez-Ortiz and J. J. de Pablo, The Journal of Chemical Physics , 014108 (2015).[75] A. H. Barnett, J. Magland, and L. af Klinteberg, SIAM Journal on Scientific Computing41