Effect of the surface charge discretization on electric double layers : A Monte Carlo simulation study

The structure of the electric double layer in contact with discrete and continuously charged planar surfaces is studied within the framework of the primitive model through Monte Carlo simulations. Three different discretization models are considered together with the case of uniform distribution. The effect of discreteness is analyzed in terms of charge density profiles. For point surface groups, a complete equivalence with the situation of uniformly distributed charge is found if profiles are exclusively analyzed as a function of the distance to the charged surface. However, some differences are observed moving parallel to the surface. Significant discrepancies with approaches that do not account for discreteness are reported if charge sites of finite size placed on the surface are considered. © 2007 American Institute of Physics. DOI: 10.1063/1.2741520


I. INTRODUCTION
Ionizable interfaces are present in many natural and industrial systems.In particular, numerous physical, chemical, and biological processes are governed by electrostatic interactions between charged colloids, surfactant monolayers, functionalized latex and oxide particles, polyelectrolytes, etc. [1][2][3] In all these systems, the appearance of electrical forces is due to the surface charges developed on the interface by different charging mechanisms.Therefore, the presence of this surface charge implicates a distribution of ions around the charged surface which is normally termed electric double layer ͑EDL͒.
For many years, the description of the EDL has been performed on the basis of the Poisson-Boltzmann ͑PB͒ equation, which becomes the cornerstone of the Gouy-Chapman model, also referred as the classical EDL theory.Therein, surfaces are assumed to be smooth and uniformly charged, whereas ions are considered to be point charges immersed in a dielectric continuum.Although the classical approach has been successfully used for the cases of model systems, these previous assumptions can give rise to erroneous predictions in the behavior of real ionizable interfaces.On the one hand, the inclusion of the ion size in the description of EDLs with multivalent ions tends to produce differences in the ionic profile which become more significant as magnitudes such as the surface charge, the ionic strength, and the ionic valence increase. 4These discrepancies explain, at least in part, several counterintuitive phenomena in nature such as the condensation of DNA, whose application in the formation of DNA complexes constitutes nowadays a feasible alternative to viral-based methods in gene delivery. 5On the other hand, since the ionizable surface charge in real systems arises from diverse mechanisms ͑such as the dissociation of chemical groups on the surface or the chemical binding of electrolyte ions͒, its nature is usually nonuniform.2][13] Other experimental systems studied as a planar ionizable interfaces are the titration of a variety of metal oxide and hydroxide particles, for example, goethite ͑␣-FeOOH͒, hematite ͑␣-Fe 2 O 3 ͒, and silica ͑see Ref. 14 and those quoted therein͒.Spherical, monodisperse polystyrene latex particles are also studied with discrete ionizable functional groups on the surface. 15Here, the functional groups are likely to be distributed randomly on the surface.At any rate, the surface heterogeneity can have considerable effects on electrokinetic properties or colloidal forces. 16,177][28][29] Computational methods can also become a valuable tool to shed light on surface heterogeneity effects.In fact, models of discretely charged surface have already been included in some of the simulations cited above.Recent examples of this ͑in the framework of the primitive model͒ are the works performed by Taboada-Serrano et al. 30 and Meyer et al. ͑among others͒. 31Taboada-Serrano et al. have employed Monte Carlo ͑MC͒ simulations to study ion distributions in the presence of a discretely charged planar surface, reporting interesting results on electrolyte mixtures, whereas Meyer et al. have analyzed the influence of the relative position of the networks of surface sites on the stability of these charged interfaces.To our knowledge, however, explicit comparisons between different models of discrete surface charge and with the case of continuous and uniform distribution are rather scarce so far.An exception concerning colloidal forces is the recently published study carried out by Khan et al. 32 They have shown that the interaction forces at short distances between surfaces characterized by a regular array of discrete charged sites could be quite different to those originated by surfaces with a uniformly smeared-out charge.In this case, the discrete nature of the surface charge is responsible for the differences in the ionic profile near the surface.
Accordingly, the main goal of this work is to deepen into the influence of the surface charge discretization on the ionic profile of a planar EDL.MC simulations have been performed to study different models of discretely charged planar surfaces in the presence of mono-and multivalent finite size ions.These models are also compared to the classical one in which a continuous distribution is assumed.In this way, the deficiencies of classical approaches ͑involving a uniformly smeared charge and point ions͒ are clearly revealed.Besides the customary ionic profiles as a function of the distance to the surface, the distribution of charge in the lateral directions ͑parallel to the plane͒ is also looked into.In this sense, the use of helpful computer resources has allowed us to increase the number of particles in the simulations as compared to that used in previous related works.
The paper is organized as follows.First, the primitive model used here and its variations are outlined.Some technical details of simulations are also provided.Then, the ionic distributions are analyzed as a function of the distance to the charged plane ͑the so-called z profiles͒.The behavior of the local charge density near the surface is also probed.Finally, some preliminary results about the diffuse potential are presented and some conclusions are highlighted.

II. MODEL AND SIMULATION METHOD
The Monte Carlo method has been used to obtain the equilibrium properties of the EDL.The calculations have been carried out in a canonical ensemble for a collection of N + positive ions and N − negative ions confined in a rectangular prism of dimensions W ϫ W ϫ L. A squared impenetrable charged wall is located at z = 0, whereas another impenetrable wall without charge is placed at z = L. Periodic boundary conditions and the minimum image convention were employed in the x and y directions. 33The simulation box contains an ionic mixture according to the concentration of the bulk electrolyte solution together with an excess of counterions neutralizing the surface charge.
Simulations were performed in the framework of the socalled restricted primitive model; thus ions were modeled as charged hard spheres of radius a and the solvent ͑water at 25 °C͒ is simply included by means of its dielectric constant.Thus, the interaction energy between two ions was calculated according to where Z i is the valence of ion i, e is the elementary charge, 0 r is the permittivity of the dielectric continuum, r ij is the relative position vector, r ij = ͉r ij ͉ is the distance between ions i and j, and d is the ion diameter ͑d =2a͒.Monovalent, divalent, and trivalent ions have been considered in the simulations.The corresponding hydrated radii a were taken as 0.36, 0.42, and 0.48 nm, respectively. 4he long-range corrections were introduced following the method employed by Boda et al. 20 based on the original idea developed by Torrie and Valleau. 18Accordingly, each ion has an associated charged plane of infinite dimensions outside the simulation box, parallel to the charged wall.A similar plane is also associated with every charged site of the charged wall.This approach does not require distribution functions obtained during previous steps and, thus, guarantees the random character of the MC simulations.
The charged wall represents a portion of an ionizable interface ͑colloidal particle, polyelectrolyte molecule, etc.͒ with a given distribution of ionized functional groups.Four models have been considered to represent such distribution ͑Fig.1͒: a continuum model ͑CONT͒ and three discrete models characterized by the radius and position of a square array of monovalent charged sites that represent the ionized functional groups.For the DISC1 model, the charged sites are points situated just on the surface.For the DISC2 model, the sites are charged hard spheres with a finite radius of 0.3 nm whose centers are also situated on the surface.The value of 0.3 nm corresponds to the estimated nonhydrated radius of three usual functional groups ͑carboxylate, sulfonate, and phosphonate͒.These estimations have been obtained from the sum of the van der Waals radius of the oxygen ͑0.15 nm͒ plus the C v O, S v O, and P v O bond distances ͑0.13, 0.14, and 0.15 nm, respectively͒. 34Bond distances have been taken from a standard force field of a molecular dynamics package. 35The values obtained for these three functional groups were 0.28, 0.29, and 0.30 nm, respectively.For the DISC3 model, the charged sites are located inside the wall at a distance of 0.3 nm from the surface.In this way, certain situations in which the charge groups are located under the interface can be modeled.Such charged groups could arise from inner defects in molecular structures.This is the case of many minerals, clays, oxides, and other compounds that can undergo isomorphous substitution, which means that structural ions are replaced with ions of valency of 1 less than the original.For example, a silicon atom ͑+4͒ in clay may be replaced by aluminum ͑+3͒, producing a surface with a net negative charge.A similar effect can be observed in the substitution of sulfate by chloride in a crystal lattice.
For these three surface models, the surface site-ion interaction energy has been calculated as in the ion-ion case.For the continuum model, however, the interaction energy of ion i with the charged wall is where r i is the position vector of particle i, z i its distance to the surface, and is the surface charge density of the charged wall.
The surface charge densities considered in this work have been −0.02,−0.04, −0.08, −0.12, −0.18, and −0.24 C m −2 .In the discrete models, they correspond to a distance between adjacent sites of 2.83, 2.00, 1.42, 1.16, 0.94, and 0.82 nm, respectively.These charge densities are similar to those found in a great variety of interfaces ͑including surfactant monolayers, latex, and oxide particles with ionizable functional groups such as carboxilate, amino, and aldehyde͒. 1It is observed that the maximum density of charged groups varies from 0.95 nm −2 for carboxilate to 1.55 nm −2 for aldehyde groups, corresponding to a mean separation between adjacent groups from 0.8 to 1.1 nm.The exact values depend on the degree of the ionization of the monolayer.For a common group density of 1 nm −2 , which corresponds to a group separation of 1 nm, and a maximum degree of ionization, the charge density is −0.16C m −2 .Smaller values are expected when the degree of ionization of the interface decreases.In any case, all of these are in perfect concordance with those employed in the present simulations.
The ionic strengths considered in this work have been 1.0, 0.6, 0.1, 0.06, and 0.01M.Accordingly, the number of ions in the simulation box varied from 1500 to 4169 particles, from 1089 to 3083 particles, and from 1344 to 2671 particles for monovalent, divalent, and trivalent counterions, respectively.The dimensions of the rectangular prism, W and L, vary from 9.9 to 28.8 and from 16.5 to 257.6 nm, respectively.The number of ions in the simulation box was fixed during any given simulation.Monodimensional z profiles were obtained from simulation runs of 3 ϫ 10 4 sampling steps per particle.For bidimensional profiles 1.2ϫ 10 6 sampling steps per particle were used for averaging.In all cases, the first 500 sampling steps per particle were employed to equilibrate the systems.Then, configurations were saved for analysis purposes every 10 sampling steps per particle.The simulations were performed by a code developed in C under a LINUX operating system in a 28 CPU cluster.
As mentioned above, our simulations have been carried out in the framework of the primitive model.In our opinion this is an acceptable choice if one is interested in looking into the EDL structure with realistic ion sizes rather than the effect of the molecular nature of the solvent.[38][39][40]

III. RESULTS AND DISCUSSION FOR Z PROFILES
The effect of discreteness is first analyzed for point surface sites.A comparison between DISC1 and CONT models is done in Sec.III A. Then, the effect of finite site on the charged surface sites is analyzed in Sec.III B comparing CONT model with DISC2 and DISC3 models.It should be noticed that, dealing with discretizated surface charge, the ionic local concentrations depend on the three coordinates, i.e., x, y, and z.However, in this section we will focus on the behavior in the direction perpendicular to the charged plane.Therefore, z profiles might be interpreted as a mean value averaged over the directions parallel to the plane.FIG. 1. ͑Color online͒ Pictures of the four models for the surface charge distribution: ͑a͒ continuous charge distribution ͑CONT model͒, ͑b͒ discrete distribution of point charge sites over the surface ͑DISC1 model͒, ͑c͒ discrete distribution of spherical charge sites of radius a over the surface ͑DISC2 model͒, and ͑d͒ discrete distribution of charge sites inside the surface ͑DISC3 model͒.

A. Point surface sites
First, we will compare results obtained for the DISC1 and CONT models beginning with the case of symmetric monovalent electrolytes.The charge distribution data will be analyzed in terms of where c i ͑z͒ is the local molar concentration of i ions and N A is Avogadro's number.In Fig. 2͑a͒, this quantity is plotted as a function of the distance to the charged surface ͑normalized by the ion diameter͒ for different salt concentrations and a rather low surface charge density ͑−0.04 C m −2 ͒.As can be seen, the results obtained for DISC1 and CONT models are practically identical.In both cases, increases in approaching the charged plane and reaches a maximum just at the distance of closest approach ͑z = a͒.This obviously means an accumulation of ͑positive͒ charge near the ͑negative͒ surface.
For ionic strengths ranging from 0.01 to 0.6M, the net charge is always positive.For 1M, however, is slightly negative around z = 1.5d.In that region, therefore, coions predominate.This situation, which cannot be predicted by the classical Poisson-Boltzmann theory, is the result of a strong accumulation of counterions just near the plane, which might eventually overcompensate its charge.
This phenomenon, known as charge inversion, can be more clearly revealed by studying the total charge accumulated up to a distance z ͑including the immobile surface charged sites͒, given by where z 0 = 0 for CONT, DISC1, and DISC2 and z 0 = −0.3nm for DISC3.The accumulated charge corresponding to data in Fig. 2͑a͒ are shown in Fig. 2͑b͒.As can be concluded, this function is clearly negative from 0.1 to 0.1M.In contrast, for 1M acc exhibits a maximum with positive values ͑i.e., charge inversion͒.This finding should be highlighted since it confirms that this phenomenon can occur for low surface charge densities.In this case, entropy plays a key role. 41According to Messina et al., the entropy of the solution ͑entropy of mixing͒ is decreased by enlarging the size of the salt ions, which enhances interparticle correlations.On the other hand, the interface provided by the macroion leads to an increase of microion density close to the macroion and promotes there certain lateral ordering, even in the absence of strong electrostatic coupling, similar to a freezing phenomenon.Figures 3͑a͒ and 3͑b͒ show the behavior of and acc , respectively, for a considerably charged surface ͑−0.24C m −2 ͒.Concerning the comparison between the DISC1 and CONT models, an almost perfect coincidence must be reported again.Regarding charge profiles, a strong accumulation of counterions near the surface is observed, but the most remarkable fact is a little pronounced maximum in for z Х 1.5d.This could be interpreted as the formation of a second layer of counterions in response to packing problems at the surface since, in a perfect close-packed structure, such a layer would be placed at ͑1/2+ ͱ 3/4͒d.In Figs. 4 and 5 the transition from low to high surface charge densities is specifically explored for two salt concentrations, 0.1 and 1.0M, respectively.DISC1 and CONT models provide coincident data once more.Apart from that, Figs.4͑a͒ and 5͑a͒ prove that the above-mentioned maximum is a feature of extremely charged systems, confirming that it is a layering effect under such conditions.In relation to this, Fig. 5͑b͒ shows how the distance required to neutralize the surface charge significantly grows for −0.24 C m −2 , which is also a consequence of this second layer of counterions.
In any case, the existence of such layering can be clearly revealed by examining the wall-ion distribution functions g͑z͒ defined as where c i * is the molar concentration of i ions at the bulk solution.g i ͑z͒ profiles are shown in Figs.4͑c͒ and 5͑c͒ for three surface charge densities, −0.04, −0.12, and −0.24 C m −2 ͑small, intermediate, and large, respectively͒.For −0.24 C m −2 , the maximum at z Х 1.5d clearly indicates that counterions undergo packing problems at 0.5d, which forces other ions to pack themselves in a second layer.The counterion layering at high surface charge density has been previously reported and extensively studied for 1:1, 1:2, and 2:1 electrolytes by Lamperski and Bhuiyan. 42What is more, these authors have tested a mean field PB approach coupled with an ionic exclusion volume term.With increasing the ionic strength up to 1.0M, the maximum corresponding to this effect is even more pronounced as the reader can infer comparing Figs.4͑c͒ and 5͑c͒.However, the latter figure also shows that, at intermediate and low surface charge densities, the ionic species located at 1.5d are the coions.This feature was not observed in the case of 0.1M ͓Fig.4͑c͔͒.The accumulation of coions at this position is a result of the charge inversion caused by the layer of counterions just near the charged surface.In the case of −0.04 C m −2 the coion distribution function shows how this species reaches the surface.As the number of counterions required for charge inversion in this case is relatively small, some coions are allowed to approach the surface.
At this point, we will examine some representative results for 2:1 electrolytes ͑and, consequently, divalent counte-rions͒.In Fig. 6, ͑z͒ and acc ͑z͒ are plotted for a series of surface charge densities ranging from −0.02 to C m −2 and a salt concentration of 0.333M, which is equivalent to an ionic strength of 1M ͑previously analyzed for 1:1 salts͒.As can be inferred, the equivalence between DISC1 and CONT models remains in this new situation.Regarding ͑z͒, the progressive evolution to a little pronounced maximum at 1.5d with increasing the surface charge density ͑discussed above for Figs. 4 and 5͒ seems to have disappeared for divalent counterions.Instead, we find again negative net charges in that region, suggesting the existence of a layer of coions, which is clearly corroborated for all surface charge densities by the maxima observed at z ϳ 1.5d in the wall-ion distribution functions in Fig. 6͑c͒.This indicates that divalent counterions do not have so many packing problems near the charged plane.Since the valence has been doubled, the local counter- ion concentration near the surface is lower compared with the monovalent case.Figure 6͑c͒ shows coion ͑instead of counterion͒ layering at z ϳ 1.5d, which was also reported by Lamperski and Bhuiyan. 42Apart from that, Fig. 6͑b͒ also illustrates how charge inversion is much more intense for divalent counterions.
In Fig. 7, snapshots for the DISC1 model of 1:1 and 2:1 electrolytes are shown for a great surface charge density ͑−0.24C m −2 ͒.It can be seen that at this surface charge den-sity, the saturation in the first layer of counterions is present for 1:1 ͓Fig.7͑a͔͒ and 2:1 ͓Fig.7͑b͔͒ electrolytes.In addition, Fig. 7͑a͒ also illustrates the above-mentioned layering effect responsible for the maximum around 1.5d for 1:1 electrolytes in Fig. 5͑a͒.Apart from that, Fig. 7͑b͒ shows how the amount of coions near the surface grows with doubling the counterion valence.
To end this subsection, some illustrative results for trivalent counterions will be considered.In Fig. 8, ͑z͒ and acc ͑z͒ are plotted for the same series of surface charge densities and a salt concentration of 0.166M, also equivalent to an ionic strength of 1M.First, it should be stressed that the predictions of DISC1 and CONT models are again practically identical.This clearly suggests that the equivalence of both models is a completely general feature for ͑averaged͒ z profiles.Regarding the net local charge, this property exhibits minima of negative values around 1.5d, whereas the accumulated charge shows a noticeable inversion at z Ϸ d, which increases with the surface charge densities ͑without exceptions for 3:1 electrolytes͒.This indicates that packing problems near the charged surface have completely disappeared for trivalent counterions, which is somehow logical since the neutralization of the plane involves a lesser number of them.Anyhow, this would allow the formation of a coion layer just behind the first counterion layer and as a consequence of charge inversion.However, for z ranging from 1.8 to 3d, just behind the region of overcompensation, acc ͑z͒ recovers its initial negative sign.This oscillatory behavior is due to alternate layers of co-and counterions.
The formation of a coion layer for 3:1 electrolytes is again revealed in Fig. 8͑c͒ by maxima at z ϳ 1.5d for the three surface charge densities studied.This figure also shows how counterions are expelled ͓g͑z͒ Ͻ 1͔ by coions from the second to a third layer ͑located at 2.7d approximately͒.This alternation of co-and counterions is responsible for the above-mentioned oscillatory behavior of acc ͑z͒.

B. Charged surface sites of finite size
Now, some representative results for models with surface sites of finite size will be discussed.First, the case for which the sites are charged hard spheres with a 0.3 nm radius located just on the surface ͑DISC2͒ will be analyzed.Figures 9͑a͒-9͑c͒ show the ͑z͒ profiles for different surface charge densities, an ionic strength of 1.0M, and 1:1, 2:1, and 3:1 electrolytes, respectively.The most remarkable feature is the significant disagreement between the predictions of DISC2 and CONT models, in contrast with the case of point ions ͑DISC1͒.For moderate and high surface charge densities, ͑z͒ exhibits a maximum shifted up to z Ϸ 0.8d due to the finite size of surface sites and mobile ions.In addition, the height of this maximum is significantly lesser, whereas its width is greater.In other words, this model of discretization markedly modifies the charge distribution near the surface.This should be taken into account in developing certain models in which the continuous approach is preferred due to its simplicity ͑e.g., adsorption models͒.What is more, it should be stressed that these discretization effects also take place for 1:1 electrolytes, for which the Poisson-Boltzmann theory is assumed to be valid.Consequently, one should be especially careful if precise values of certain parameters ͑e.g., binding constants͒ were obtained from experimental data by a fitting procedure based on this classical approach. 43,44Such theory might lead to wrong estimations.Apart from that, it must be pointed out that the discrepancies due to surface charge discretization also involve some disagreement between the predictions of these models ͑DISC2 and CONT͒ for acc ͑z͒ ͑not shown in these figures͒.The snapshots of Fig. 7 also show differences between DISC1 and DISC2 models, for example, in the reduced number of coions present in the same region near the surface.The comparison between 1:1 and 2:1 electrolytes through snapshots for DISC2 ͓Figs.7͑b͒ and 7͑c͔͒ leads to conclusions similar to those commented above for the DISC1 case.
Finally, the situation in which the surface sites of finite size are located at z = −0.3nm ͑DISC3͒ will be briefly commented.The results obtained then ͑see also Fig. 9͒ considerably resemble those reported for point surface sites before.In particular, the agreement with the continuous model is again outstanding.This is reasonable since, being the charged sites inside the surface, the hard-sphere interactions between them and mobile ions are not operative.

IV. CHARGE DENSITY IN THE PROXIMITY OF SURFACE SITES
The preceding section has been mainly devoted to the behavior of the charge density in the direction perpendicular to the surface ͑z͒.As we have seen, only one of the surface discretization models yields significant differences with the case of continuous charge distribution.In this section, however, the local net charge density, ͑x , y , z͒, will be examined in the vicinity of the charged surface sites as a function of x and y for z = a and z =2a.
Figure 10 shows the comparison of the charge density of the solution around a surface charge site from DISC1, DISC2, and DISC3 models for a 1:1 electrolyte, 0.1M in the presence of a surface charge density of −0.12 C m −2 .These data have been obtained by averaging over all sites of the surface.In all cases, the charge site is located in the center of the x and y axes.The cross sections corresponding to diagonal directions ͑e.g., y = x͒ are also plotted.It can be seen that the ion distribution differs greatly in the proximity of the surface depending on the model considered.For DISC1, a maximum of charge density is observed around the position of the site at z = a.For z =2a, nevertheless, this maximum has almost disappeared.For the DISC2 model, at z = a ͓Fig.10͑b͒, left͔, there is no density around the position of the charge site because of the site-ion volume exclusion.However, maxima are observed in the position that corresponds to the central point of the line connecting surface sites ͑that is, the corners of the corresponding cell͒.For z =2a, some ordering along x and y directions remains.However, its structure is different.The map together with its cross section reveals that the charge is now preferentially located around the surface sites ͑instead of the cell corners͒ but less structured, since the difference in height between the peaks and the cen- tral depression is considerably smaller than for z = a ͑the reader can compare cross sections͒.Finally, for the DISC3 model, nearly plane profiles are obtained at z = a and z =2a.For this model, charge sites are located inside the surfaces.Thus, the magnitude of the electrostatic interactions of these sites with the solution is smaller than for the other models.Consequently, charge density distributions along x and y directions are practically independent of the position of the surface charge site, at least for 1:1 electrolytes.
Figure 11 shows ͑x , y , z͒ as a function of x and y for z = a and z =2a but now for a 2:1 electrolyte.For z = a, conclusions similar to those pointed out for Fig. 10 can be reported.However, ordering effects along x and y directions are now stronger, as expected.Particularly, in the DISC2 model, counterions are again located at the corners of each site cell forming a two-dimensional square lattice ͑as surface sites form͒.The excluded volume effect might be partly responsible for the two-dimensional crystal mentioned above.It should be stressed, however, that electrostatic forces ͑more intense for divalent ions͒ might also contribute to the formation of this lattice.In fact, multivalent counterions can form two-dimensional crystals even in the absence of surface charge discretization as a result of strong electrostatic correlations. 45,46or z =2a, the slight depression previously reported for DISC2 and monovalent salts ͓see Fig. 10͑b͔͒ has, however, disappeared for 2:1 electrolytes.In order to justify this behavior, the reader should keep in mind that hydrated divalent ions are larger than the monovalent ones.Consequently, we are examining ͑x , y , z͒ further from the charged surface than in the monovalent case and the ordering effects along x and y directions tend to vanish with increasing z.Finally, it should be noticed that even though the z profiles are nearly identical for the CONT, DISC1, and DISC3 models, the differences found in the charge distribution in the lateral directions might involve a noteworthy effect in the behavior of electrostatic forces at short distances between particles.In particular, it could be important to take them into account in the correction of binding constants through the surface potential, due to the fact that the mean field approximation commonly used in the determination of surface potential cannot be appropriate for these cases. 44

V. DIFFUSE POTENTIAL: SOME PRELIMINARY RESULTS
Having simulated the ionic distributions, the calculation of other properties such as the diffuse potential should be addressed.In the case of a uniform surface charge density, this task can be easily carried out from the z profiles as follows: 4,20 For discretized surface charge distributions, however, ion densities also depend on x and y.Therefore, this way is not available.One could then apply other standard methods such as those based on the definition of electrostatic potential.
Unfortunately, for the charged plane, such methods present serious difficulties mainly related with the fact that the origin of the electrostatic potential for ions and the charged sheets employed to include long-range corrections in the energy are not the same.At any rate, we would like to give a preliminary idea of the effects of surface charge discretization on the diffuse potential.In order to do that, the dependence on x and y could be neglected and a first estimation of the diffuse potential would be done from Eq. ͑6͒.The reader should, however, be aware of the approximations involved in this procedure.
Here, d was actually computed by applying a procedure equivalent to Eq. ͑6͒, but averaging the resulting diffuse potential for different Monte Carlo ionic configurations N conf and for different chosen values N z max of the maximun distance z max , from which it is reasonable to suppose bulk conditions, Six different values of z max around the half-length of the simulation box were selected.For each value of z max chosen, 2500 MC ionic configurations are used to calculate an estimation of the diffuse potential.The integral term in Eq. ͑6͒ is computed as a summation extended over all k ions, N ions , with z k ഛ z max,j , being z max,j the maximun distance chosen in each case.W 2 is the surface area considered in the simulation box.This procedure only considers the z distance of each ion to the surface and neglects the dependence on x and y ͑as mentioned before͒.In this sense, the obtained value could be seen as an approximate mean diffuse potential.Figure 12 illustrates the behavior of the diffuse potential as a function of the surface charge density for DISC1 and DISC2 and for an ionic strength of 1.0M of 1:1, 2:1, and 3:1 electrolytes.As can be easily inferred, the model of surface charge discretization has considerable effects on this property, particularly at moderate and high surface charge densities.For DISC2, the potential curves are shifted toward the region of negative values.This is logical since this model considers finite surface sites ͑unlike DISC1͒.Consequently, the neutralization of surface charge by adjacent counterions would be lesser ͑as compared to DISC1 or CONT͒.

VI. CONLUSIONS
In this work, several models of discretely charged surfaces have been compared by computing the ionic distributions of the EDL in the presence of these charged surfaces.Moving perpendicularly to the charged plane and comparing with the case of continuous distribution, no differences are reported if point surface sites ͑DISC1͒ are considered.In other words, we have shown that the discretization of the surface charge as point sites on its own does not modify the z profiles.In this sense, the differences in the local charge density due to such discretization are compensated in averaging in x and y directions.For finite size sites located just on the plane ͑DISC2͒, however, significant discrepancies appear even for 1:1 electrolytes, for which the classical EDL approach ͑involving a continuous surface charge͒ is assumed to be applicable.As a result of the excluded volume effects, the ionic profiles near the surface, ͑z͒, are not so abrupt and the ion concentrations can be noticeably smaller ͑for moderate and high charge densities͒.
Besides the z profiles, the local net charge density near the surface sites has also been studied for all the discretization models.This sheds light on the singular behavior of the DISC2 model.In particular, we have seen that the predictions of the three discretization models do differ.For DISC1 and DISC3 models, counterions tend to locate over the surface sites.On the contrary, for DISC2 these ions situate be- tween surface sites.In any case, these configurations must be considered to understand the influence of surface charge discretization on colloidal stability.

FIG. 2 .
FIG. 2. ͑Color online͒ Charge density of a 1:1 electrolyte solution ͑a͒ and its integration ͑b͒ as a function of the distance to a −0.04 C m −2 charged surface for CONT ͑lines͒ and DISC1 ͑symbols͒ models.Profiles for five different electrolyte concentrations are shown.

FIG. 4 .
FIG. 4. ͑Color online͒ Charge density of a 1:1 electrolyte solution of 0.1M of ionic strength ͑a͒, its integration ͑b͒, and wall-ion distribution functions ͑c͒ as a function of the distance to a charged surface for CONT ͑lines͒ and DISC1 ͑symbols͒ models.Profiles for several different surface charges are shown.In ͑c͒ solid symbols stand for counterions whereas open symbols indicate coions.

FIG. 5 .
FIG. 5. ͑Color online͒ Dependence of the charge density of a 1:1 electrolyte solution of 1.0M ionic strength ͑a͒, its integration ͑b͒, and wall-ion distribution functions ͑c͒ as a function of the distance to a charged surface for CONT ͑lines͒ and DISC1 ͑symbols͒ models.Profiles for several different surface charges are shown.In ͑c͒ solid symbols stand for counterions whereas open symbols indicate coions.

FIG. 10 .
FIG.10.͑Color online͒ ͑x , y͒ charge density obtained from DISC1 ͑a͒, DISC2 ͑b͒, and DISC3 ͑c͒ models for a 1:1 electrolyte solution of 0.1M ionic strength at distances z = a ͑left column͒ and z =2a ͑right column͒ from a −0.12 C m −2 charged surface.Each plot shows a squared portion of the ͑x , y͒ function centered around a surface charged site.The monodimensional charge density profiles, ͑r͒, corresponding to a diagonal section of ͑x , y͒ surfaces are also shown ͑d͒.