A Multi-scale Monte Carlo Method for Electrolytes

Artifacts arise in the simulations of electrolytes using periodic boundary conditions (PBC). We show the origin of these artifacts are the periodic image charges and the constraint of charge neutrality inside the simulation box, both of which are unphysical from the view point of real systems. To cure these problems, we introduce a multi-scale Monte Carlo method, where ions inside a spherical cavity are simulated explicitly, whilst ions outside are treated implicitly using continuum theory. Using the method of Debye charging, we explicitly derive the effective interactions between ions inside the cavity, arising due to the fluctuations of ions outside. We find that these effective interactions consist of two types: 1) a constant cavity potential due to the asymmetry of the electrolyte, and 2) a reaction potential that depends on the positions of all ions inside. Combining the Grand Canonical Monte Carlo (GCMC) with a recently developed fast algorithm based of image charge method, we perform a multi-scale Monte Carlo simulation of symmetric electrolytes, and compare it with other simulation methods, including PBC+GCMC method, as well as large scale Monte Carlo simulation. We demonstrate that our multi-scale MC method is capable of capturing the correct physics of a large system using a small scale simulation.


I. INTRODUCTION
Monte Carlo (MC) simulation of Coulomb many body systems is a difficult problem [1,2], mainly due to the long range nature of Coulomb interaction. In the most naïve simulation strategy, one would confine a collection of ions inside a box with hard walls, and compute the total energy by adding up all pairwise interactions. This leads to two obvious difficulties: 1) The artifacts of hard wall propagate into the bulk of the system, with the characteristic scale set by the Debye length. This renders large amount of simulation data useless. The situation becomes worse in the dilute limit. 2) The computational complexity of each Monte Carlo cycle is of order of N 2 , with N the system size. This makes simulation of large size systems practically impossible. Obviously, both difficulties have their roots in the long range nature of Coulomb interaction.
Adoption of periodic boundary conditions (PBC) restores the translational symmetry and hence eliminates all the boundary effects. Furthermore, one can take advantage of periodicity and compute the long range part of Coulomb energy in Fourier space, using the method of Ewald summation [3]. (For a pedagogical introduction, see [2].) This reduces the computational complexity from N 2 to N 3/2 , which is still not fast enough for large scale simulations. For molecular dynamics (MD) simulations, one can use particle-mesh method and fast Fourier transform [4,5] to reduce to complexity further down to N log N . Unfortunately, these techniques are not applicable for Monte Carlo simulations. This is probably the main reason why MD has become the more popular method for charged systems, even when one is only insterested in static properties.
By periodically replicating the system, one introduces an infinite array of image charges for each charge belonging to the system under study. These images are unphysical and break the rotational symmetry. Furthermore, the total charge inside the simulation box must vanish, for otherwise the summation over images would diverge. This leads to an unphysical constraint of charge neutrality inside a finite volume. In Fig. 1(a), we show the average charge density around a fixed ion, obtained using different simulation methods. Whilst large scale MC simulation (STD) gives a clean form of screened Coulomb, small scale simulation (with system size three Debye lengths) using PBC (PBC+GCMC) gives substantially different results. The deviation, which grows with the distance to the test ion, is caused by the unphysical image charges. Similarly, internal energies calculated using PBC also deviate substantially from the standard results, shown in Fig. 1(b). These artifacts are quantitatively less important for larger systems, and become invisible for infinite systems. In principle, one can always simulate sufficiently large systems so that these artifacts become insignificant, or correct these artifacts for every microscopic configurations during the simulation. This however generically leads to waste of computational resources, or slow-down of simulation processes.
As an alternative to PBC, one may use a linear continuum theory to model the influence of the subsystem outside the simulation domain, whilst particles inside the domain are simulated explicitly. For a dipolar system, the subsystem outside responds to the dipoles inside the cavity and exerts a field on the latter, which are conventionally called the reaction field. FIG. 2: Schematic illustration of the multi-scale reaction potential grand canonical Monte Carlo. If this simulation method work well, the boundary of simulation cavity must behaves as a virtual one, i.e., the simulated system must behave as ions inside a virtual spherical domain in an infinitely large system. field boundary condition is named after Onsager [6], who first used this method to calculate the dielectric constant of dipolar fluids. Similar multi-scale ideas, i.e. treating dipolar molecules in near and far fields using different methods, however can be dated back to Clausius [7] and Mossotti [8]. For a discussion, see the monographs on dielectrics by Fröhlich [9], and by Battcher [10]. Born [11] also used a similar idea to calculate the solvation (free) energy of ions in solvent. We note that in the study of ionic systems, the term reaction potential would be more appropriate than reaction field, since it is the electrostatic potential, rather than field, that couples to the charges. The same method can be generalized to molecules with arbitrary charge distributions. Simulation methods using multi-scale strategy have been studied by many authors [12][13][14][15][16][17][18][19][20]. These methods were designed for study of either dipolar systems or ionic systems. The simulation cavity can be either fixed and common for all simulated particles, or can be moving and individual for each particle. The subsystem outside the cavity is always assumed to obey Poisson equation (for dipolar systems) or the linearized Poisson-Boltzmann equation (for ionic systems). Finally, it appears that in most of the previous works, the reaction-field modeling is used in molecular dynamics simulations (MD), even though in principle, it is applicable for Monte Carlo simulations as well.
To implement the multi-scale simulation, one must solve the continuum theory for each simulation step. The classical Kirkwood expansion [21], though straightforward, is not efficient because of its slow convergence. The image charge method, which was first discovered by Neumann [22] for polar systems, and have been generalized to ionic systems [23][24][25], transforms the Kirkwood expansion into a line integral, and therefore can be efficiently computed. Recently, we improved the image method [26] (using Mellin transform technique) so that it is applicable for Poisson-Boltzmann theory.
Regardless of the rich history, there are still a few con-tributions that we can make to improve this multi-scale simulation strategy. Firstly, in order to mimic an infinite system using a finite and simulation domain (which is the essential task of almost all computer simulations), the total charge inside the simulation cavity must be allowed to fluctuate. This means that the numbers of positive and negative ions must be allowed to fluctuate independently. A "grand canonical Ewald" simulation where the ions are inserted and deleted in neutral pairs can not capture the physics of charge fluctuations. In this work, we shall discuss in detail the artifacts due to suppression of charge fluctuations. Secondly, all previous works assume without proof that the subsystem outside the cavity is described by the linearized Poisson-Boltzmann theory (LPB) for ionic systems. We shall derive the correct linearized continuum theory and show that it reduces to the naïve LPB only in the dilute limit. For non-dilute electrolytes, the continuum theory may be considerably complicated. Furthermore, for asymmetric electrolytes, there is also a constant cavity potential that acts on all ions inside the cavity. Thirdly, we shall combine the multi-scale GCMC with a recently developed image charge method [26] to simulate symmetric electrolytes. We find that multi-scale simulation of a small system (with linear size of only three Debye length) is capable of capturing the physics of an infinite system. For example, as shown in Fig. 1(a) and Fig. 1(b), both the screening charge density around a fixed ion and the average internal energy per ion computed using our multi-scale GCMC (RP+ GCMC) agrees remarkably well with the standard results (STD). By contrast, to achieve similar precision in simulations using PBC, the system would have to be at least as large as ten times Debye length. The remaining of this paper is organized as follows. In Sec. II, we analytically discuss LPB with PBC, and demonstrate two sources of artifacts of PBC, i.e., periodic image charges, and constraint of charge neutrality inside each unit cell. In Sec. III, we start from the basic principle of statistical mechanics, and systematically derive the multi-scale MC simulation strategy. This derivation made it clear that one should always use grand canonical ensemble for the subsystem inside the cavity, and that, for asymmetric electrolytes, there is a constant cavity potential acting on every ions inside. In Sec. IV we outline the numerical implementation of the multi-scale strategy using grand canonical MC simulation, as well as that of other simulation methods that are used for comparative studies. In Sec. V, we simulate the symmetric primitive model, and compare the simulation results using different methods. We demonstrate how our multi-scale method captures correctly the correlation effects between ions, whilst the Ewald summation methods fail to do so. Finally in Sec. VI we draw the conclusion remarks.
To simulate large systems, one need to use fast algorithms based on multipole expansion, such as the octtree algorithm [27], to expedite the computation of electrostatic energy. A Monte Carlo simulation strategy that combine multi-scale modeling, grand canonical en-semble, image charge techniques, and oct-tree algorithm speeding-up can achieve a computational complexity of order of N log N for one Monte Carlo cycle, and is free of artifacts due to PBC. It therefore constitutes a competitive alternative to molecular dynamics simulation methods for electrolyte systems. This will be discussed in a future work.

II. ARTIFACTS DUE TO PERIODIC BOUNDARY CONDITIONS
Consider N positive ions q and N negative ions −q inside the simulation cube, which is centered at the origin and has length L, and volume V = L 3 . To simplify the analysis, we shall only focus on symmetric electrolytes in the dilute regime, where linearized Poisson-Boltzmann theory is applicable. We shall for the moment assume that all ions that point-like, and consider the case of finite ion sizes towards the end of this section. We periodically replicate the system, so that both the charge distribution and the electrostatic potential become periodic functions.
We fix a positive ion q at the origin, and let all other ions fluctuate according to the equilibrium Gibbs measure. We shall calculate two quantities: 1) the screening charge density ρ q (r) around the test ion, and 2) the correlation potential ψ, i.e., the mean potential acting on the test ion due to all other ions. These two quantities are calculated in the classical Debye and Hückel theory (DH) [30] for an infinite electrolyte. What we are doing here is to generalize DH to a finite system with PBC. The total internal energy of the system is related to ψ via Now according to the linearized Poisson-Boltzmann theory (LPB), the average number densities of mobile positive and negative ions are related to the mean potential φ(r) via the Boltzmann factor: where we have linearized the exponentials in the second steps. For a physical system with infinite size, one usually requires that the mean potential φ(r) vanishes in the infinity. Then according to Eqs. (2),c ± are the number densities of positive/negative ions in the infinity, i.e., the bulk ion densities [33]. This argument is no longer applicable if periodic boundary condition is imposed, since the mean potential φ(r) does not vanish as infinity. For periodic systems, the most convenient choice is that the mean potential φ(r) vanishes identically when integrated over the unit cell: This is equivalent to the "conductor boundary condition" that are popularly used in Ewald summation method, which requires that the Fourier transform of electrostatic potential does not have a k = 0 component. Now if we integrate Eqs. (2) over the unit cell with volume V , the LHS are the total numbers of mobile positive/negative ions in the unit cell, which are N −1, and N , respectively, whereas the RHS arec ± V respectively. This allows us to determine the parametersc ± as Let us now invoke the Poisson equation, use Eqs. (2), and expand to the first order in φ. We obtain where l = (L x n x , L y n y , L z n z ) are lattice vectors, n x , n y , n z are integers, andκ is the inverse Debye length: The first term in RHS of Eq. (5) is due to the test ion and its periodic images, whilst the second term corresponds to a uniform and negative charge density, arising due to the fact there is one more negative mobile ion than positive mobile ions, see Eq. (4). This in turns originates from the artificial constraint of charge neutrality inside each primitive cell. Integrating both sides over a unit cell, one easily see Eq.
(3) is indeed satisfied. Note that inverse Debye lengthκ in Eq. (6) is also slightly different from its usual form in free space. Eq. (5) can be easily solved subject to PBC: which can be easily shown to satisfy the condition Eq. (3).
Taking the Laplacian of this equation, we find the total average charge density: where the first term is due to the test ion and its images, whilst the second term is sum of the screening charge densities around these ions. Eq. (8) (averaged over all orientations of r) is plotted in Fig. 1(a) as PBC+DH, which agrees remarkably well with MC simulations with PBC (PBC+GCMC). The correlation potential ψ acting on the fixed ion can be easily obtained by subtracting off the part due to the test ion itself from Eq. (7), and take the local limit r → 0: The first term is just the correlation potential in free space as predicted by Debye-Hückel theory, and the other two terms constitute the artifacts due to periodic boundary conditions, both of which vanish as the system size L goes to infinity (with the Debye length fixed).
If the ions are not point-like, but are hard spheres with diameter d (the primitive model), the above results need to be modified properly. The resulting theory is considerably more complicated. If the ion density is not very high, however, we may consider corrections due to hardcores as small perturbations, same as the artifacts due to PBC. Then all we have to do is to correct the first term in Eq. (9), i.e., replace it by the correlation potential of a hard sphere ion, which was worked out by Levin and Fisher some time ago [28,29]. This leads to This result is plotted in Fig. 1(b) (PBC+DH), which agrees reasonably well with the MC simulations with PBC if the ion density is low. The agreement becomes increasingly worse as the density increases, indicating that the approximation underlying Eq. (10) becomes increasingly inaccurate.
In Fig. 1, we have chosen the size of simulation box to be approximately three Debye length, in order to highlight the artifacts due to PBC. In many simulations, the size of simulation box is chosen to be much larger than the Debye length. In this case, the second term in Eq. (10) is exponentially small comparing with the third term. The artifacts is therefore completely due to the last term. The relative error (in the dilute regimeκa → 0) due to PBC is the given by To make this relative error less than 1%, one would have to choose the size of simulation box to be larger than ten times Debye length. This may be very difficult to achieve in many simulations. By strong contrast, using the multi-scale GCMC simulation method, we can faithfully simulate an infinite system using a simulation cavity with radius only three times Debye length.

A. Effective Interaction
We shall study the primitive model of electrolytes, where the solvent (water) is modeled implicitly as a dielectric medium with relative permittivity / 0 = 80, whilst ions are modeled as hard spheres with a point charge at the center. Furthermore, the permittivity of the hard spheres is assumed to be the same as that of the solvent. We introduce a spherical simulation cavity, which is already schematically illustrated in Fig. 2. The subsystem A inside the cavity consists of N charges {q 1 , . . . , q N } with positions denoted by X = {x 1 , . . . , x N }, whereas the subsystem B outside the cavity consists of M charges The total Hamiltonian of the overall system consists of three parts: where H A and H B are the Hamiltonian of isolated subsystems A and B: The first term in each bracket is the long-range Coulomb interaction and the second term w(x − y) is the short range interaction, which is assumed to be independent of the ion species . For the primitive model, w(x − y) is just the pairwise hardcore repulsion: The last part in the RHS of Eq. (12), H AB , represents the interaction between subsystems A and B: We note that the short range part of H AB in Eq. (15) becomes nonzero only if the distance between x i and y α is smaller than d, which can happen only if both particles are very close to the cavity boundary, hence it has no influence on the ions that are at least one ion-diameter away from the boundary. We shall therefore ignore this short range interaction between A and B. Accordingly, we shall exclude a thin spherical shell near the cavity boundary for collection of simulation data. H AB can then be written into the following form: where ϕ(x i , Y) is the potential acting on the charge q i at x i due to all ions in B, in the specific micro-configuration Y = {y 1 , . . . , y M }: The grand canonical partition function of the overall system can be written as where the traces Tr A and Tr B represent integrations of all the variables as well as summation over particle numbers: Here N s and M s are the numbers of mobile particles of s-th species in subsystems A an B respectively, whereas µ s are their chemical potentials, and Λ s are some microscopic length scales. We integrate out the variables in subsystem B, and obtain an effective theory for the subsystem A: where δH A represents the additional effective interactions of subsystem A as mediated by subsystem B: Therefore we can study the subsystem A only with an effective Hamiltonian and the influences of the subsystem B on A are automatically taken care of. We will find that δH A can be approximately calculated and the result assumes a simple and physically transparent form. We emphasize that Eqs. (21a) and (18) are grand canonical partition functions. Hence the ion numbers for each species in A and B are all stochastic variables. Now it is well known that for large systems (such as subsystem B), all ensembles are equivalent to each other. But for small systems (such as A), different ensembles are inequivalent. It is therefore important to use grand canonical ensemble, rather than canonical ensemble for A in our multi-scale modeling. We shall demonstrate this point by comparing simulations results using GCMC and CMC in Sec. V.

B. Method of Debye Charging
To calculate the effective interaction δH A , we shall use the Debye charging method [30]. Let us scale charges of all ions q i in subsystem A by a factor λ: so that Eq. (16) becomes Eq. (21b) then becomes: Taking derivative of Eq. (25) with respect to λ, we find where the λ-dependent average O λ is defined as The physical significance of ϕ(x i , Y) λ is therefore the average potential at x i (inside the cavity, occupied by an ion λq i ) due to all ions in B. Note that the average is over the statistical fluctuations of all ions outside, with the locations of all ions λq i inside fixed. To calculate this average, we treat λ as a small parameter and use the linear response theory. That is, we expand Eq. (27) to the first order of λ: Note that the average O 0 in RHS is defined by Eq. (27) with λ set to zero, i.e., with the interaction between two subsystems switched off: Let us define a kernel χ(r, r ) in terms of the connected correlation function of ϕ(r, Y) as follows: Setting O = ϕ(x i , Y) in Eq. (28) and using Eqs. (29), (16), we obtain: The first term in the RHS is the average potential at x i in the absence of any charges in A, generated by an isolated subsystem B. We shall call this term the cavity potential Φ cav (r), following the terminology of Onsager [34]. Evidently it satisfies the Laplace equation inside the cavity, and is invariant under arbitrary rotation. It then follows that Φ cav (r) must be a constant in the whole cavity: Furthermore, for symmetric electrolytes, the cavity potential vanishes identically due to charge inversion symmetry. For asymmetric electrolytes, Φ cav is generally non-vanishing. This term has been missed by all previous works on reaction-field modeling of charged systems. The second term in the RHS of Eq. (30) is linear in λ. It therefore arises due to the linear reaction of subsystem B to sources charges {λ q 1 , . . . , λ q N } in A. We shall therefore call it the reaction potential [35].
Finally, combining Eqs. (31) and (30), and substituting them back into Eq. (26), integrating the latter over the charging parameter λ from 0 to 1, we obtain the effective interaction of A mediated by B: The correlation function χ(r, r ) is intimately related to the electrostatic Green's function of the subsystem B. To see this, let us insert a test ion with charge q at r inside the otherwise empty cavity, and calculate the average total potential at r, which may be inside or outside the cavity. This potential is the linear superposition of a part due to the source charge q, and another part due to all ions in subsystem B. The later has the form given by Eq. (30), to the first order in q (with λ set to unity, of course). Hence the total mean potential at r is given by The sum inside the bracket describes the linear response of mean potential at r due to a unit point charge at r , and therefore is the electrostatic Green's function: Substituting this back into Eqs. (22) and (32), we find the effective Hamiltonian H A eff : where the first term is the effective interaction between different ions, and the second term is the self-energy of ions mediated by subsystem B.

C. Linearized Poisson-Boltmann Theory
The Green's function G(r, r ) encodes the linear response properties of the subsystem B, that is an electrolyte with a spherical void (the empty cavity). More precisely, it is the incremental mean potential at the field point r due to the insertion of a unit test charge at the source point r . If r is inside the cavity, G(r, r ) obviously satisfies the Poisson equation: If r is outside the cavity, it satisfies a linear integrodifferential equation: (36b) The kernel α(r, r ) is generally unknown. At the level of first order perturbation, however, we may approximate α(r, r ) by the kernel that corresponds to a uniform electrolyte, which only depends on the difference r − r , because of the translational symmetry. The Fourier transform of this kernel was calculated analytically recently, both for symmetric and asymmetric electrolytes [32]. The results are still quite complicated.
In this work, we shall only study symmetric electrolytes in the dilute regime, where LPB is applicable in the bulk. This entails two essential simplifications for the effective Hamiltonian Eq. (35). Firstly the cavity potential vanishes identically. Secondly, the kernel α(r, r ) has the simple form κ 2 δ(r − r ), where κ is the inverse (bare) Debye length, given by with c the average ion density. Consequently, Eq. (36b) reduces to the linearized Poisson-Boltzmann equation: − ∆ r G(r, r ) + κ 2 G(r, r ) = 0, r ∈ B. (36d) We shall deal with the cases of dense electrolytes and asymmetric electrolytes in a future work. We also need to determine the boundary conditions satisfied by the Green's function. At infinity, it clearly satisfies the free boundary condition: On the cavity interface, r = R, G(r, r ) and its normal derivative are continuous: where r → R ± mean that the field point r approaches the interface from outside/inside, respectively.
Our method therefore works as follows. Firstly we find the Green's function by solving Eqs. (36a), (36d), subject to boundary conditions (37), then use Eq. (34) to find the reaction potential χ(r, r ), and finally use Eq. (35) to find the effective Hamiltonian.

IV. NUMERICAL IMPLEMENTATION
In this section, we discuss the numerical implementation of our multi-scale MC method, as well as other simulation methods that we use for comparison.

A. Computation of the Green's Function
In one of our previous works [26], we discussed an efficient image charge method for solving Eqs. (36a) and (36d), subjected to boundary conditions Eqs.(37). Here we briefly summarize the results in order to make this work self-contained. Let r and r be, respectively, the source point and the field point. We obtain the reaction potential χ(r, r ) in the form of Kirkwood series [21]: where θ is the angle between r and r , P n (x) the Legendre polynomials, R the cavity radius, u = κR, and M n (u) = (n + 1)k n (u) + uk n (u) nk n (u) − uk n (u) , with k n (u) the modified spherical Hankel functions: Eq. (38) has a useful integral representation: where the vector t is parallel to r and can be written as Its magnitude t runs from r K = R 2 /r to infinity. Eq. (41) has an intuitive physical interpretation: it is the potential at r due to a fictitious straight-line of image charges which starts at the Kelvin point r K = R 2 r /r 2 , and extends all the way to infinity. As shown in reference [26], the line charge density λ(t) [36] is related to the Mellin transform of M n (u) (as a function of n), and is an oscillatory function of t. To further reduce the computational cost, we truncate the line integral and discretize the line image using Gaussian quadrature. χ(r, r ) can then be expressed by M point-like image charges and a few correction terms: where M image point charges are located at r m , with magnitude q m . The correction terms are due to the truncation of the line integral. Numerical tests showed that 4 images and L = 1 corrections provide result with error less than 1% in computing the self-energy of a charge. All relevant details of q m , r m and χ l can be found in reference [26]. For not too small r, Eq. (42) is much more efficient than the Kirkwood expansion, Eq. (38). If the source point approach the center of the sphere, r R, the Kirkwood expansion converges rapidly. We can use leading term of the Kirkwood expansion to approximate the Green's function: . (43)

B. Grand canonical Monte Carlo
The system parameters used in our multi-scale simulation are as follows. The radius of simulation cavity is fixed to be 100Å. The total ion number (including both species ) varies from 15 to 40. The corresponding range of Debye length is between 24 − 50Å. Hence the system size is not much larger than the Debye length, and it is a nontrivial matter to cancel the influences of boundary and restore the physics of an infinite system. All ions have hardcore diameter d = 7.5Å.
We simulate the subsystem A with an effective Hamiltonian Eq. (35), with the reaction potential computed using Eq. (42). As emphasized in the preceding section, we must use grand canonical Monte Carlo to simulate this system. Use of canonical Monte Carlo would suppress fluctuations of total charge inside the cavity and therefore leads to substantial errors. We will demonstrate this point in great detail below.
The probability density function of a microstate with N ions {q 1 , . . . , q N } at X = {x 1 , . . . , x N } is given by the grand canonical Gibbs distribution: (44) Note that the grand canonical partition function Ξ is dimensionless, whilst the dimension of p(X) is L −3N . To perform numerical simulation, however, it is mandatory to deal with discrete probabilities. Therefore we discretize the simulation domain, with elementary length unit |δx|. The probability of a discrete microstate is then The Grand Canonical Monte Carlo simulation consists of three steps: displacement, insertion and removal [1]: (1) Displacement. We select an ion randomly, and let its charge and position be q i and x i . The part of total energy of the system that depends on q i is given by Then we choose a destination x i for q i according to a flat probability density function defined in a cubic cell with a given size and centered at x i , and let the new energy be E i . The change of the total energy of the system is The displacement is accepted with probability Evidently, if any hardcore constraint is violated in the destination state, ∆E = ∞, and the displacement is rejected with probability one. If the center of an ion across the boundary of simulation domain in an attempted move, it is also rejected. In another word, the boundary behaves as a hard wall to the center-of-mass of all ions. It is easy to see that this transfer probability satisfies the detailed balance with respect to the equilibrium Gibbs distribution Eq. (45): π(X) acc (X → X ) = π(X ) acc (X → X) . (49) (2) Insertion.
To mimic the fluctuations of particle numbers in an open system, we need to insert and remove particles. In a micro state with N s particles of species s, we first randomly select a position x in the cavity according to a uniform probability distribution. The latter can be generated by introducing three independent random variables η 1 , η 2 , η 3 that are uniformly distributed in the interval [0, 1], and express the spherical coordinates of the insertion point x as It is straightforward to verify that the Jacobian ∂(x, y, z)/∂(η 1 , η 2 , η 3 ) is a constant independent of η 1 , η 2 , η 3 , where x, y, z are the Cartesian coordinates. We now choose a species s randomly and insert an ion q s at x with probability acc (N s → N s + 1) = min 1, where E i is the energy change due to the insertion operation, given by Eq. (46), whilst µ * s are the effective chemical potential of the sth species , given by (3) Removal.
In a micro state with N s particles of species s, we randomly choose a particle with (a randomly chosen) species s, and remove it with probability where −E i is the energy change of the removal operation, with E i again given by Eq. (46).
It is straightforward to verify that the probabilities of insertion and removal satisfy detailed balance with respect to the equilibrium grand canonical distribution Gibbs distribution Eq. (45): The factor |δx| 3 /V in the LHS is the probability of choosing one particular point inside the (already discretized) simulation domain, whereas the factor 1/(N s + 1) in the RHS is the probability of choosing a particular particle with species s.

C. Comparison with Other Simuation Methods
Three other simulation methods are used to compare with our multi-scale reaction potential GCMC simulation method (RP+GCMC). To demonstrate the artifacts due to periodic boundary conditions, we use the popular Ewald summation GCMC method (PBC+GCMC) to conduct a small scale simulation. To demonstrate the artifacts due to suppression of charge fluctuations, we apply a small scale simulation on reaction potential model using canonical Monte Carlo (RP+CMC). Finally, we conduct a large scale simulations with system sizes at least ten Debye lengths to obtain standard results (STD), with which all other simulations are compared. To make the comparison meaningful, we adjust parameters such that the ion densities in the bulk are equal in all simulations. Furthermore, all small scale simulations have equal size of simulation domains.
(1) Grand canonical Ewald summation MC with Periodic Boundary Conditions (PBC+GCMC). We choose cubic simulation domain with volume equal to 4πR 3 /3. To avoid divergence when summing over periodic image charges, the total charge in the simulation box must vanish in every micro state. Furthermore, one still need to specify some "boundary conditions" at infinity. We shall use the popular "conductor boundary condition", where the average potential inside the simulation box vanishes, and there is no dipolar term in the total free energy. For details, see the textbook by Frenkel and Smit [2]. This choice is consistent with the condition Eq. (3) that we use in the analysis of LPB with PBC.
The total ion number however can fluctuate according to the standard grand canonical ensemble theory. In accordance with these constraints, ions are inserted and deleted at random as pairs [1].
(2) Multi-scale Reaction Potential Canonical MC simulation (RP + CMC). This simulation is similar to RP+GCMC, with the only difference that the ion number of each species is kept constant.
(3) Large Scale Canonical MC simulation with Hard Wall Boundary Conditions (STD). We conduct a large scale canonical simulation using a spherical simulation domain with radius at least ten times of Debye length, and with hard wall boundary condition. To void the artifacts due to boundary effects, only ions that are more than five Debye lengths away from the boundary are used for data collection. The simulated system typically contains about 2500 to 5000 ions. To speed up the computation of total electrostatic energy, we use the recently developed octtree algorithm [27]. We calculate various physical quantities via this method and use them as standard results for all later comparisons.

V. RESULTS AND DISCUSSIONS
We compute various physical quantities using three small scale simulations, and compare them with the standard results (STD). Disagreement with STD indicates the simulation strategy is incapable of capturing the correct physics of infinite systems.

(1) Average ion densities inside cavity
We measure the average ion density in the simulation cavity as a function of the radial distance, shown in Fig. 4. If the multi-scale modeling is faithful, the resulting average ion density must be independent of radius. In another word, non-uniformity of ion density indicates that the artifacts of cavity boundary have not been properly cancelled. In Fig. 4, we see that RP+GCMC simulation yields a much flatter density profile than RP+CMC. For the latter case, there is a systematic tendency that the ion density is higher near the center than near the boundary, with an overall variation of about one percent. This is caused by the artificial constraint of charge neutrality inside the simulation cavity.
Note that both RP+GCMC and RP+CMC show a thin boundary layer (with thickness comparable with ion diameter) where the ion densities vary rather abruptly. This is due to the short scale depletion effects of the hard wall boundary, and shows up actually in all simulations with hard wall boundary conditions. Even though correction of these depletion effects are straightforward, it is more convenient just to exclude the boundary layer completely for data collection.
This comparison clearly demonstrates that we should use grand canonical, instead of canonical, ensemble when implementing the multi-scale simulation strategy.
(2) Correlation potential We measure the mean potential acting on an ion due to all other ions, which is usually called the correlation potential. Because we are in the low density regime, the Debye-Hückel theory is applicable. Indeed, as shown in Fig. 5, both RP+GCMC and STD agree well with prediction of DH. By contrast, RP+CMC shows substantial deviation from the other results. This is again due to artificial constraint of charge neutrality inside the cavity.
This correlation potential is not measurable in PBC+GCMC simulation, because the summation over potential due to all periodic images is divergent.

(3) Screening charge density around a test ion
We fix an ion inside the simulation cavity and measure the screening charge density around, as a function of the distance to the fixed ion. This charge density can be obtained using data of pair correlation functions inside the whole simulation cavity. The results have already been presented in Fig 1(a). Again, RP+GCMC agrees with STD up to high precision, whilst both PBC+GCMC and RP+CMC yield results that are substantially different.
Theoretically, LPB predicts that the charge density has the form of screened Coulomb: Hence if we plot rρ q (r) v.s. r in log scale, we should obtain a straight-line, with slope given by the inverse Debye length. As shown in Fig 1(a), this agrees remarkably well with the simulation results using RP+GCMC. The charge densities obtained using both RP+GCMC and STD can be used to determine the screening length as a function of ion density. The results are shown in Fig. 6. Also shown there is the prediction of screening length by LPB, Eq. (36c), which agrees rather well with simulation results. This shows that LPB is indeed applicable for the ion densities studied in this work.

(4) Internal energy per ion
We use different methods to compute the internal energy of the system per ion. The results are shown in Fig 1(b). The total internal energy of the system can be easily calculated in terms of the correlation potential acting on each ion: E = 1 2 i q i ψ i . Again our multi-scale GCMC simulation yields the same results as the large scale simulation. Both multi-scale CMC and PBC+GCMC yield incorrect results for the internal energy. In the same figure, the solid line (DH) is the prediction of Debye-Hückel theory in open space (the first term in RHS of Eq. (10)), which agrees with STD and RP+GCMC very well. The dashed line is the prediction of Debye-Hückel theory with PBC as well as leading order correction of hardcore taken into account (Eq. (10)). It appears that Eq. (10) takes into account most of the artifacts in PBC+GCMC simulation if the density is sufficiently low. As the density increases, Eq. (10) deviates from PBC+GCMC more significantly.

(5) Total charge inside the simulation cavity
As a stringent test of the faithfulness of our multi-scale method, we fix a positive ion at a distance r from the origin, and measure the total charge inside the cavity. The geometry is illustrated in Fig. 3. The total mean charge inside the cavity is nonzero, because screening can not be perfect in a finite volume. In the framework of LPB, the total charge density is given by Eq. (55). Hence the total charge inside the cavity can be easily calculated: This result, shown as the solid solid curve in Fig. V, agrees remarkably well with STD and RP+GCMC. Again, this indicates that 1) LPB works well in the density regime under study; 2) our RP+GCMC simulation method faithfully captures the physics of an infinite electrolyte. By contrast, for RP+GCMC, this total charge is identically zero, due to the condition of charge neutrality.
(6) Ion densities near a charged colloid Finally to test the applicability of our multi-scale simulation strategy for inhomogeneous systems, we insert a uniformly charged colloid with charge Q 0 = 20e and radius R 0 = 10Å at the center of the simulation cavity, and measure the local potentialφ(r, q) acting on an ion fixed at r > R 0 , due to both the charged colloid as well as all other ions. We set the dielectric constant of the colloid to be the same as that of the solvent, so that there is no image charge effect. Note thatφ(r, q) also depends on the ion diameter d. We do not show this dependence in order to avoid cluttered notations.
Note that the potentialφ(r, q) is different from the average potential φ(r) in the absence of the test ion. We shall call the former the conditional mean potential, and the latter the unconditional mean potential. Both quantities can be directly measured using simulation data. The difference between them arises because the test ion influences the distribution of other ions. The difference betweenφ(r, q) and φ(r) encodes essential information about ionic correlations in this inhomogeneous system. We shall numerically computeφ(r, q) and clarify its relation with φ(r). This demonstrates that our multi-scale reaction potential GCMC is capable of capturing the ionic correlations faithfully in an inhomogeneous system.
The conditional mean potentialφ(r, q) can be expanded in terms of the charge q: φ(r, q) = φ 0 (r) + q φ 1 (r) + O(q 2 ). (57) It is important to realize that φ 0 (r) is not the same as the unconditional mean potential φ(r). Their difference is caused by the hardcore of the test ion. Nevertheless, in the far field regime where LPB is applicable, we can show that φ 0 (r) and φ(r) are proportional to each other, and both assume the form of screened Coulomb potential: φ 0 (r), φ(r) ∝ e −κr r .
Hence φ 0 (r), φ 1 (r) can be computed directly using simulation data. In Fig. 8(a) and 8(b), we plot r φ 0 (r) and r [φ 1 (r) − φ 1 (∞)] as a function of r in log scale. The far field asymptotics of both functions indeed behave as screened Coulomb, with their characteristic decay lengths differ by a factor of two. This demonstrates that in our multi-scale GCMC simulation, the influences of the boundary of simulation cavity are properly cancelled, so that the simulation data can be used to analyze subtle ionic correlations. Also the fact that log (r φ 1 (r)) scales linearly with r in the far field indicates that the slow function χ(r) is indeed changing very slowly. The potential of mean force (PMF) of a test ion (with positive charge q) can be obtained from the local potential Eq. (57) via the method of Debye charging: U (r, q) = q 0 φ (r, q) −φ(∞, q) dq + U (r, 0).
The constant of integration U (r, 0) is the free energy cost of bringing a neutral particle from the bulk to r. In a symmetric electrolyte, U (r, 0) decays twice as fast as φ 0 (r), and therefore makes no contribution to the leading order far field asymptotics of U (r, q). Now the average density ρ + (r) of positive ions is related to the PMF U (r, q) via the following well known Gibbs distribution. Hence the total charge density is ρ(r) = qρ 0 e −βU (r,q) − e βU (r,−q) ≈ −2q 2 β φ 0 (r). (63) That is, the deviation of ion density from its bulk value also behaves as a screened Coulomb in the far field. In Fig. 8(c) we plot the corresponding simulation results and show that indeed ρ(r) obeys the law of screened Coulomb, with the same characteristic decay length 1/κ.

VI. CONCLUSION AND ACKNOWLEDGEMENT
We have demonstrated the artifacts due to periodic boundary condition in the simulation of electrolytes. These artifacts are caused by the periodic image charges and the constraint of charge neutrality inside simulation domain, both unphysical from the perspective of mimicking an infinite system. One can always avoid these artifacts by simulating a much larger system and throwing away data from the region near the boundary. This leads to waste of computational resource, but is frequently done in simulations. Our multi-scale reaction potential GCMC provides a more efficient alternative. Combining reaction potential modeling and grand canonical Monte Carlo, this method cancel most of the electrostatic arti-facts introduced by the boundary of simulation data, and hence faithfully capture the physics of an infinite system even in a very small scale simulation (with linear size about three Debye lengths). If we combine this multiscale GCMC method with appropriate fast algorithm for evaluation of electrostatic energy (such as the oct-tree algorithm), we arrive at a method for large scale MC simulation of electrolytes, that will constitute a competitive alternative to the current methods based on molecular dynamics. This will be explored in a future publication.
We thank NSFC (Grants No. 11174196 and 91130012) for financial support, the High Performance Computation Centre (Π) at Shanghai Jiao Tong University for computation resources. We also thank Wei Cai for interesting discussions.