Comparison of simulation methods for charged systems of slab geometry

Using the specific model of a system of like charged ions confined between two planar like charged surfaces, we compare the predictions for the energy and density profile of four simulation methods available to treat the long range Coulomb interactions in systems periodic in two directions but bound in the third one. Monte Carlo simulations demonstrate unambiguously complete agreement between the results obtained with these methods where the potential between charges is solution of Poisson's equation in the simulation cell with adequate boundary conditions. The practical advantages of the different methods is assessed.

The standard way to treat long range interactions (Coulomb, dipolar, Yukawa near the Coulomb limit) in simulations of bulk systems is to replicate the basic simulation cell periodically in all three directions of space and apply an Ewald summation technique (EW3D) [1,2]. However, in many situations of electrochemical, biological or technological interest, as, for instance, electrolyte solutions between charged surfaces, charged lipid bilayers in water, suspensions of colloids between glass plates, clays [3][4][5], magnetic thin films [6], Wigner crystals [7] etc.., the system is finite in one direction thus necessitating a modification of the convential Ewald method. For systems periodic in two dimensions but bounded along the third one an Ewald summation method (EW2D) for electrostatic interactions has first been given by Parry [8,9] and was later rederived by various authors [10][11][12][13][14]. Unfortunately, the practical use of the EW2D sum is hampered by the occurrence, in the reciprocal space term, of a double sum over different particles which, due to the complicated way the bounded coordinates enter the expression, cannot be reduced, as for the EW3D case, to a sum of order N, a circumstance which renders the method computationally expensive. Not surprisingly, only few calculations using EW2D have been reported to date mainly to test the validity of more approximate approaches [15,16]. These calculations involve tabulation and inteerpolation of the potentials on a three-dimensional grid.
Various "time saving" schemes have been proposed to bypass the computational burden of the EW2D method. The purpose of this paper is to compare some of these methods for the specific case of a system of charged particles confined between two charged planar surfaces separated by a distance of the order of several particle diameters. The methods that have been chosen are those which satisfy exactly the laws of electrostatics (perfect screening, Green's theorem etc..) and do not present numerical instabilities difficult to control. Thus the charged sheet method proposed by Torrie and Valleau [17,18] and its modification by Boda et al [19] have been discarded. Indeed the distribution of the point charges of the ions outside the simulation cell, located in the image cells of the latter, is represented approximately by a set of uniformely charged planar sheets. An other method, proposed by Lekner [20,21] is based on rewriting the sum of forces acting on an ion by a second ion and its periodic images as an absolutely converging infinite sum. The energy is obtained by integration of the force. A shortcoming of the method is the need, to estimate the sum with a given precision, to retain a number of terms depending strongly on the relative position of pairs of particles. Because of this technical reason, emphasized in more detail in ref [22], the method was not included in this investigation.
The four methods we have considered, the first three of which have already appeared in the literature, are EW3D [15,16,23,24], a method developed by Hautman and Klein [25] (HK), the method of hyperspheres (HSG) [26] and that of concentric spheres (CS) described below.
A way to treat the long range Coulomb interactions is to place the slab at the centre of a parallelepipedic simulation box which has dimension perpendicular to the slab much larger than the width of the latter and apply the EW3D method. In this way one effectively simulates an infinite number of parallel slabs. Provided the region of empty space separating the slabs is sufficiently large one expects the influence of the periodic images on the behaviour of the system to become negligible.
If the distance between the surfaces confining the particles is small compared to the extension of the system along the directions parallel to the surfaces (narrow slab), the lateral distance s between pairs of particles will be much larger than the distance z normal to the slab surfaces and it will be appropriate to expand the Coulomb pair interaction in powers of z/s and apply an Ewald sum to the in-plane component of the interaction. Such an approach was followed by Hautman and Klein (HK) [25].
An attractive way to investigate the behaviour of Coulomb particles between charged surfaces which circumvents the cumbersome Ewald sums is to use a closed space, e. g. a hypersphere in four-dimensional Euclidean space [26][27][28]. Not only is the electrostatic po-tential, solution of Poisson's equation on the hypersphere obtainable in simple closed form [27], thus avoiding approximation of the interaction potential (as for instance in EW3D due to truncation of the direct and reciprocal space sums), also the number of operations necessary to calculate the distances between particles is reduced with respect to Euclidean space with concomitant speed up of the simulation procedure. This geometry has been applied previously to the study of the attraction between two like charged surfaces neutralized by solvated counterions in an endeavour to comprehend the stability of charged lamellar materials such as clays and cement [27,29] and to the study of the effective interaction of charged colloidal particles confined between like charged plates [30].
A system of charged particles occupying the region between the outer and inner surfaces of two concentric spheres is also a suitable arrangement to describe, in the limit of sufficiently large radii of the spheres, the system of charged particles confined between two planar surfaces. This geometry has the advantage that the interaction between charges and between charges and surfaces is the usual Coulomb potential; however it may require use of a large number of particles to render the effect of curvature of the confining surfaces small.
The purpose of this paper is to check, by means of Monte Carlo simulation, the agreement between the results of the four aforementioned methods for the energy and density profile of the system described above consisting of N like-charged ions between parallel like-charged surfaces separated by a distance h. The ions, modelled as hard spheres of diameter d bear a charge q at their centre and each surface, of area S, a uniform charge density σ.
Such a system can be viewed as a crude model for lamellar liquid crystals formed by ionic amphiphiles [31] or charged lamellar materials like clays or cement [27]. The convergence of the results to their thermodynamic limit is estimated by performing simulations with an increasing number of particles keeping the distance between surfaces fixed. The speed of this convergence is obviously an important criterion for appraisal of the relative practical interest of the four methods.
Expressions for the energy calculated with the different methods will be given in Sect.
II together with details of the Monte Carlo simulations. Results for the energy and density profiles obtained with the four methods are compared in Sect. III and conclusions drawn in Sect. IV.

A. EW3D
In this method a square slab of ions of thickness h is placed at the centre of a simulation box having dimension L z normal to the surfaces much larger then the lateral dimension L and the system is extended periodically in the three directions of space. The slab surfaces are located at z = ±h/2 perpendicular to the z-axis. The closest approaches of the ions to the impenetrable surfaces are therefore z = ± 1 2 (h − d). To evaluate the total energy of the system which is the sum of the contributions from the ion-ion U ii , ion-surface U is and surface-surface U ss interactions it will be convenient to associate to the ions a uniform neutralizing background (filling the whole simulation cell) of density −Nq/V and to the charged walls a uniform background of density −2σS/V = −2σ/L z where V = L z L 2 = L z S is the volume of the simulation cell. Because of electroneutrality of the system Nq + 2σS = 0, the backgrounds cancel out exactly and will not contribute to the total energy.
The contribution of the ions and their background to the total energy is given by In these equations r ij = r j − r i where r i and r j are the positions of particles i and j.
The term − 1 2 C w is the self-energy given by In Eq. (1) ν is a vector of components (Ln x ,Ln y ,L z n z ) (n x , n y , n z integers) and the prime in the sum over ν means that the terms i = j must be omitted when ν = 0. The wave-vectors G which enter the reciprocal space contributions to the energy have components 2πn x /L, 2πn y /L and 2πn z /L z In our calculations the sum on lattice vectors extends over all G subject to |n x | ≤ 6, |n y | ≤ 6, |n z | ≤ 12 and |n| ≤ n max = 12. The α parameter which governs the rate of convergence of the real-and reciprocal-space contributions was taken sufficiently large so that only the terms with ν = 0 had to be retained in Eqs. (1) and (2) The electrostatic energy between the ions and the two charged walls including their compensating backgrounds is given by where V s (z) is the electric potential associated with the electric field E s generated by the two charged surfaces and theit backgrounds. Due to the symmetry of the simulation cell, this electric field is normal to the charged surfaces and depends only on the zcoordinate. Furthermore, as a consequence of the periodic boundary conditions, it must and, by integration, for the potential, choosing V s (0) = 0 (any additional constant would In the periodic system the potential created by the charged surfaces is thus parabolic. The potential V s (z) can also be obtained by a direct integration of the Ewald potential over the two charged planes resulting in a one-dimensional Fourier series which can be summed explicitly to yield Eq. 5. Integration of V s over the volume of the simulation box leads to Finally, the interaction between the two charged surfaces (including their background) is B. Hautmann-Klein method.
In this method the slab, having the characteristics described above, has periodic boundary conditions only in the x-and y-directions. The origin of coordinates being at the centre of the parallelepipedic cell, the energy of the system is given by where r ij = r j − r i , r i ,r j coordinates of the particles i and j, and the vectors r α and r β are the position vectors of points on the surfaces S α (α = 1, 2) having constant z coordinate The sum over ν runs over all vectors of components (n x L, n y L) (n x , n y integers) perpendicular to the z-direction. The three sums in the expression for U correspond to the interactions between the particles, the interaction between the particles and the surfaces S 1 and S 2 and to that between the two surfaces S 1 and S 2 , respectively. To each of these sums, due to the infinite number of replicated cells, is associated a divergent contribution. However, if electroneutrality is taken into account these divergent terms will cancel.
The added and subtracted term represents the first M + 1 terms in the binomial expansion of 1/r in powers of z/s where s is the component of r in the plane of the surface. By introducing a convergence function h n (s; α) for each term 1/s 2n+1 the energy can be divided in a short range part and a long range part which can be evaluated in reciprocal space. In these equations s ij,ν = |s j − s i + ν|. With the choice h 0 (s; α) = erf(αs) and h n (s; α) s 2n+1 = 1 a n (2n)!
as proposed by Hautman and Klein the short and long ranged parts of the electrostatic energy take the form [25] where only the ν = 0 term has been retained (see below) and with G a two-dimensional vector of components 2π(n x /L, n y /L). In our simulations we kept terms up to M = 3 which allowed the short range part U s ii to be limited to its ν = 0 term even for relatively large values of h. Expressions for h n (s; α) up to M=3 are The attractive feature of the method is that, by writing the z 2n ij explicitly as polynomials in z i and z j , U l ii becomes a sum of terms each of which involvs products of two functions having general form In the evaluation of U l ii sum on the pairs of particles is replaced by sums on particles which, obviously, makes the computation faster. The contributions to the energy of the interactions between the ions and the surfaces and between the surfaces are easily obtained using the method of de Leeuw and Perram [11] and the identity The divergent terms of these two contributions having been eliminated through use of the electroneutrality condition, the ion-surface and surface-surface interactions contribute to the energy U by a constant term equal to 2πσ 2 V . The expression for the total energy is therefore given by

C. Hyperspherical geometry
In this method the Monte Carlo simulations are performed on a hypersphere S 3 in fourdimensional Euclidian space. Two surfaces of angular colatitudes θ N and θ S = π − θ N separated by a distance h = R(π − 2θ N ) are located symmetrically on opposite sides of the equator (see Fig. 3 of ref. [27]). N neutralizing ions of charge q are confined between the two surfaces which bear each a charge density σ. For given h and σ the aperture θ N is obtained by solving the equation Similar to the EW3D case described above, neutralizing backgrounds are associated to both the ions and the charged surfaces. As a detailed derivation of the different contributions to the potentiel energy, ion-ion U ii , ion-surface U is and surface-surface U ss can be found in ref. [27] only the final expressions are given here Here the angle α is equal to the ion radius d/2 divided by the hypersphere radius R, θ ij is the angular separation between particles i and j and q s = σS (S area of each surface).

D. Concentric spheres
In the last method we have explored, the ions occupied the region of volume V between two concentric spheres of radii r l and r m = h + r l . The external surface of the inner sphere, of area S l = 4πr 2 l , and the internal surface of the outer sphere, of area S m = 4πr 2 m , bear a charge density σ. The two surfaces being impenetrable the distances of closest approach of an ion to the surfaces are r l + d/2 and r m − d/2. The electroneutrality conditions reads The total energy of the system is readily obtained as where the three first terms in the r.h.s. represent the ion-ion energy, the energy of the ions with the charged surface S l and with the surface S m which creates in the volume V a constant potential. The three last terms correspond to the energies associated with the surface charges of S l and S m .

III. RESULTS
In our comparisons we did not consider realistic values of q and σ as for instance envisaged in ref. [27,29]. Our aim being methodological, we chose values of q and σ allowing for a compelling test of the efficiency of the four methods to predict reliably the properties of the confined charged system. In the following charges q and σ will be expressed by using reduced units of charge and distance (βq 2 /d) 1/2 and d, respectively. In these units the surface charge σ has been fixed to -1 in all our simulations and the value of q varied from 0.5 to 5. The distance between the surfaces limiting the system has been taken to be h = 5; however, in order to test the limit of validity of the HK method for large h a few simulations were The simulations were carried out in the canonical ensemble. Of the order of 10 5 MC steps per particle were performed to calculate the energy and the density profiles when N ≤ 3000.
In the CS geometry, for N ≈ 15000, 10 4 trial configurations per particle were generated to calculate these quantities. The statistical error on the energies is of the order of ±0.02 in units of kT . For the two planar geometries, EW3D and HK, the density profiles are estimated with an error of ∼ 0.5%; the statistical error, except for systematic error due to curvature effects, is of the same order for the HSG method. In the case of the CS geometry, for N ≈ 15000, the error on the energies is of the order of ±0.03 and on the density profiles of the order of 2 − 3%.
In the calculations using the EW3D method, the simulation volume has a square section of side L and an elongation along the z-axis of L z varying between 60 and 90. For h = 5 and L between 15 and 30, the influence of the value of L z on the simulation results turned out to be always negligible. Also, within this geometry, if the system confined between the surfaces has a net dipole moment, a correction to the energy proportional to the square of the dipole moment normal to the surfaces can be taken into account to remove the interaction between net dipoles in the periodically repeated slabs [16]. As the dipole moment in our system was small such a correction was not considered.
A systematic comparison of the energies for 0.75 ≤ q ≤ 5 is made in Table 1. It shows without ambiguity the convergence of the results of the four methods. It appears that the thermodynamic limit for U is obtained for q ≥ 1 as soon as the number of particles is of the order of 500 for EW3D and HK and of the order of 2000-3000 for the HSG and CS method.
However, at q ≈ 0.75, i. e. at densities larger than 0.4 at fixed h and σ, at least 10000 particles are necessary for the energy calculated with CS to agree within 1% with the other methods. This result is obviously related to the important curvature effects expected for radii r l and r m ≤ 10. Figure 2 shows, for the four methods, the variation of U when the ionic charge q varies from 0.5 to 5. A maximum is observed for q ≈ 0.8. This maximum occurs as a consequence of the lowering of q and saturation of the particles in contact with the confining surfaces due to steric effects. Indeed layers of particles parallel to the surfaces form at high density as evidenced by Fig. 3 which presents the density profiles for the different values of q. For large q, layers of highly localized particles are in contact with the surfaces. When the density increases (i. e. q decreases) layers separated by a distance of order d appear in the volume. The profiles shown in Fig. 3 are obtained by the method of HK and are indistinguishable from those given by EW3D. at q = 0.75 is in agreement with the slow convergence, in the latter geometry, of the energy to its thermodynamic limit when q < 1.
Extrapolations towards z = ±2 give the values ρ c of the profiles when the particles are in contact with the surfaces. It has been shown that for planar surfaces separated by a distance much larger than the diameter d of the particles the pressure is given by [32,33] where E c , the value of the electric field at the surface, equals 4πσ for planar surfaces. For the confined systems considered in this work this expression is not strictly valid anymore.
Nonetheless, considering it as an approximation and using the values of ρ c given in Table 2, an estimate of P can be obtained. The value of P is small and positive for q > 1 and rises to a value of 5 at q = 0.5 where hard core interactions are important. The nearly zero value of the pressure for q > 1 is clearly in accord with the absence of particles in the central part of the simulation cell (cf. Fig. 2) and complete screening of the charged surfaces by the particles in their vicinity implying an almost vanishing electric field in this same region of simulation volume. 14 and 15 for the energy, M is chosen to be 3.

IV. CONCLUSION
The convergence of the results obtained with the four simulation methods considered is clearly the main conclusion of this study. It demonstrates the strict equivalence of the various possible routes to take into account the long range of the Coulomb interactions : use of periodic boundary conditions, a closed system without boundary or simply a large system with boundaries. The equivalence is made realized only by the use of potentials preserving the laws of electrostatics and thus solution of Poisson's equation associated with the simulation cell and its boundary conditions. One can remark that the system under investigation is particularly challenging to establish this equivalence because particles all bear a charge of the same sign and screening effects therefore result only from interactions between the particles and the surfaces.
Depending on the type of system which is simulated, the advantage of using one method or the other may differ. This work clearly shows that the EW3D and HK methods apply more favourably to the study of systems confined by planar surfaces, the thermodynamic limit being reached for N ∼ 500. However, with regard to computing time, the method of hypersphere is the most efficient despite the necessity to use systems with larger values of N to make curvature effects negligible. Obviously, the CS method is unfavourable to study planar interfaces as curvature effects get small only for N > 10000 − 20000. It has been used here mainly to establish without ambiguity the equivalence between the Ewald and hypersphere "Coulomb" potentials and the usual Coulomb potential.
The calculations carried out with the CS method nevertheless show that the ionic density profiles are affected in the vicinity of a charged interface by the curvature of the latter.
Such modifications of the density profiles near planar charged interfaces may be present in solutions of inverted micelles or suspensions of charged colloids of small diameter.
The determination of the relative efficiency of the simulation methods for the calculation of the properties of confined systems will be complete only when pressure, free energy and surface tension will have been evaluated. Calculation of these quantities is in progress.