Electrical double layers (DLs) formed by ions near a charged surface are everywhere, from everyday technologies (e.g., near electrodes in batteries) to biological cells (e.g., near cell membranes and around proteins and DNA) to lab devices and materials (e.g., inside porous Nafion and nanofluidic devices), to name just a few. And while DLs have been studied for more than a century starting with the classical theory of Gouy-Chapman-Stern [1, 2, 3], DLs have become more important than ever because of new technological applications that exploit the physics and energetics of DLs. For example, the capacitive energy stored in DLs forms the basis of supercapacitors (a.k.a. electrochemical capacitors)  and similarly capacitively holding ions near a charged surface is used for desalination . Temporary protection of inner metal surfaces of industrial equipments from corrosion (e.g., in oil refinery operations) is highly affected by the properties of DLs. In these applications, the ideal ions are those that provide a high capacitance (i.e., a lot of charge stored with little applied voltage). Other potential applications like optimizing sensors that read the binding and unbinding of charged aqueous ligands  utilize ions that provide low capacitance . DLs can also be used to efficiently convert pressure into voltage in nanofluidic devices  when their dimensions become comparable to the Debye length [9, 10, 11].
This growth in the potential technological applications of DLs has also spurred a growth in theories and computational techniques to describe them [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]. The need for new theories comes from the fact that the GCS theory treats ions as dilute point charges that interact through a mean electrostatic potential, an assumption that breaks down when ion size and non-mean field electrostatics produce correlations between ions [50, 28, 51]. These correlations can manifest themselves as oscillations in the density profile of ions near electrode and in charge inversion where a layer of co-ions forms behind the first layer of counterions and can change the sign of the electrostatic potential near the electrode. Such phenomena have been unraveled through computer simulations [52, 53, 54, 50, 55, 56, 57, 18, 19, 28, 29, 36, 45]. As our new technologies progress, they are encountering these correlations and have started to exploit them .
One problem, however, has been that there is no systematic data available against which to test new theories, something that is important as some well-established theories of DLs are qualitatively incorrect . Here, our goal is to address this gap by providing high quality Monte Carlo (MC) simulation data for the primitive model of ions over a wide range of ion sizes, valences, and concentrations. The ion and electrostatic potential profiles we provide for download in Supporting Information can be used directly to better understand ion correlations and capacitance in general as well as in special laboratory and industrial procedures like electrophoresis . This data can also be used to check the range of validity of new theories of DLs so that creators and users of these theories can have confidence in when they are and are not correct.
Moreover, this data explicitly provides the capacitance of the DL in the primitive model for the cases we study since we provide the entire (discretized) electrostatic potential profile; the voltage holding the charge per unit area in the DL (whose total is negative of the surface charge because all simulations have well-established baths in the middle of the simulation cell) is the potential at the wall.
The primitive model is the simplest model of ions to include size by modeling them as charged hard spheres in a background dielectric that approximates the solvent. While all-atom simulations with explicit waters using molecular dynamics would be the ideal (see the review of Spohr  and references therein), computing power and quality of the force fields have not reached the point where one can do a systematic study like the one here. That being said, the primitive model has been used for more than 40 years , and it has provided significant insight into DLs and reproduced many experimental results, including those with ion correlations like charge inversion [63, 64, 65].
Since the seminal papers of Torrie and Valleau [52, 53, 54, 50, 55] who developed the Grand Canonical Monte Carlo (GCMC) technique for the DL geometry, MC simulations of the primitive model DL have been the gold standard because of the relatively few assumptions and approximations they use. While the ion density profiles have statistical fluctuations due to finite sampling, we have tried to minimize this by running long simulations. Our intent is to provide a comprehensive, high-quality database for as many cases as can realistically be dealt with GCMC or canonical simulations.
Model and method
The electrolyte is modeled as a mixture of charged hard spheres immersed in a solvent represented as a dielectric continuum of dielectric constant . The pair potential describing the interaction between two ions is then
where and are the valence and radius of ionic species , is the elementary charge, is the permittivity of vacuum, and is the distance between two ions. Ion diameters are denoted by .
We used the GCMC method  to simulate this primitive electrolyte model (a misleading term; this model is simple, but far from being “primitive”) confined between two charged hard walls at and along the -dimension. The interaction potential between such a charged hard wall and an ion is
where is the distance of the ion from the th surface. The surface charges, , are equal in magnitude and opposite in sign on the left and right walls. In our simulations, the left wall is always the negative one.
The cross section of our finite simulation cell in the dimensions is a square. Periodic boundary conditions are applied in these dimensions. The effect of the periodic images of the ions in the central cell are taken into account with a modified version  of the charged sheet method of Torrie and Valleau .
The GCMC method includes random insertion/deletion of neutral ion clusters ( cation and anions, where and are the stoichiometric coefficients). For this, we need to know the mean chemical potential of the salt, . This was determined with the Adaptive GCMC method of Malasics et al. [67, 68] for a prescribed salt concentration denoted by . In our study, we deal with 1:1, 2:1, and 3:1 electrolytes, so the salt concentration is equal to the cation concentration. The temperature was 298.15 K throughout this study.
The raw data obtained from the simulations are the density profiles, , computed from the average number of ions in layers of width divided by the volume of the layer (). For convenience, we report our results for the density profiles in unit of M (mol/dm), so we call them concentration profiles.
The electric field and the potential were computed from solving Poisson’s equation . The electric field was obtained by integrating the charge profile, . The integration constant was obtained from assuming that the electric field is zero outside the solution domain. The electrostatic potential was obtained by integrating once more. The integration constant was chosen so that . Then, we spatially averaged the potential over the central part of the simulation cell (around ) and shifted the profile with this value. In the simulations, we took care that a bulk region formed in the middle of the simulation cell, where and are close to constant (this is not always easy, at low concentrations, for example). The result is the potential with respect to its value in the bulk.
The width of the bin was Å in every simulation. This value fit an integer number times into the ionic radii used in this work, which makes the integration of the charge profiles more straightforward. This small value, however, results in extremely noisy profiles for low concentrations (see Fig 1). These profiles could be smoothed (by averaging into larger bins, for example, as shown in Fig 1), but we report the raw data because they contain all the information provided by the simulations. We trust the user to smooth the profiles. Despite this noise due to the fine grid, the electric field and potential profiles are smooth (see inset of Fig 1).
The simulations are long enough (200 millions sampled configurations, namely, attempted insertion/deletions or ion displacements) so potential profiles do not change considerably after running longer simulations.
When GCMC neutral particle group insertion/deletions fail to establish a bulk concentration (their acceptance ratio can be far below 0.1% in the most difficult cases), we used canonical simulations. Their advantage is that we just displace the ions individually within the simulation box, and the acceptance ratio of that kind of MC step is much higher. The disadvantage is that we need to iteratively recalculate the dimensions of the simulation cell to achieve the desired bulk concentration. The GCMC process is more convenient (if works), while the canonical iterative process is cumbersome, but sometimes it is the only way (e.g., 3:1 electrolyte with M).
We performed simulations for an extended range of state variables. The valence and diameter of the anion were fixed at the values and Å, respectively. We changed the cation attributes by varying its valence (, 2, and 3) and diameter (, 3, 6, and 9 Å).
We varied the bulk concentration in steps by orders of magnitude (, 0.001, 0.01, 0.1, and 1 M) so we cover a wide range, from extremely dilute to quite concentrated. The surface charge is varied in finer steps up to C/m (0, 0.02, 0.04, 0.06, 0.08 and 0.1 C/m). This makes it possible to compute capacitance at zero electrode charge (data not provided). We also tested quite high surface charges up to C/m (0.175, 0.25, 0.375, and 0.5 C/m), where it was possible.
Simulations were straightforward mostly, but in several cases it was difficult to establish a bulk region in the middle of the simulation cell. Therefore, we do not report results for the following cases:
M and M for the 3:1 case
M for the 2:1 case
M for the 1:1 case at high surface charges
These cases were difficult for the following reasons:
Low concentrations with their large Debye lengths are complicated because the simulation cell must be long enough ( Å) in order to fit the wide DL into it. Large cross sections, , must also be used in order to include enough ions into the cross section so we can sample ion correlations not only in the , but also in the dimensions. Sampling of possible ion configurations in this huge space is difficult when we have an otherwise dilute system.
It is especially difficult in the case of 2:1 and 3:1 electrolytes, where we insert/delete 3 or 4 ions instead of just 2. The acceptance ratio of the insertions/deletions, therefore, is very small especially at high concentrations. In most cases, however, we were able to collect enough sample from insertions/deletions except for 3:1 electrolytes at M using large cations ( and 9 Å). In these cases, we used canonical simulations.
The other especially difficult case is when the surface charge is large because the system is divided into a very high density (for the ions) region near the electrode and a low density region in the bulk.
Fig 1 illustrates the problematic case associated with low and high .
The other problematic case is large concentration ( M) with multivalent cations (especially ). In this case, the fine noise associated with low density is replaced by a coarser noise associated with high density and low acceptance ratios of particle moves (Fig 2). In these cases, a few percent of error in the reproduced bulk densities can occur. These errors have little effect on the potential profiles.
Also, a long range ordering can be observed (oscillations in profiles that extend far from the wall) that is sometimes difficult to distinguish from noise (Fig 2). The structure of the DL near the electrode is, however, well established; the potential profiles, therefore, can be calculated accurately. We can state that, in general, the potential is dominated by the region near the electrode, while the tail of the diffuse layer reaching into the bulk has a smaller effect.
The potential at the electrode is shown in Fig 3 as it is a quantity of basic interest. The potential at the cation radius (Fig 4) is closely related to the otherwise poorly-defined zeta potential that is supposed to be the potential at the slip plane that separates immobile ions attached to the electrode from those that are mobile in the diffuse layer. This value, therefore, has practical importance. A few conclusions on the basis of these figures follows:
For the 3:1 systems, the electrode potential is independent of concentrations for negative surface charges because the trivalent cations adsorbed closely to the electrode dominate the DL. The diffuse layer has very small effect.
In the 3:1 case, we see negative differential capacitance for Å (Fig 3).
An interesting maximum appears in Fig 4 for Å, which is the result of the balance between opposing effects of large ionic charge (increasing field at the ion’s surface) and large ionic size (decreasing field at the ion’s surface). Fig 2 shows the effect of these competing asymmetries in the concentration and potential profiles.
We have produced systematic simulation results for the classic model of the DL over a wide range of parameters. Some of these parameters (small concentration, large electrode charge, large ionic charge, large ionic size) have historically posed a challenge for reasons we discussed. Now, however, there is sufficient computing power and patience to overcome many of these challenges. There is also a compelling scientific need. By making this data freely available, research groups can use it directly in their work or to test new theories.
Supporting database. A compressed file that is a zipped subdirectory system organized into a tree structure. This subdirectory system contains all the concentration, electric field, and electrostatic potential profiles along with basic simulation parameters for each run. Nomenclature for directory and file names is given in the S1 File.. Available from figshare.com
DB and MV gratefully acknowledge the financial support of the National Research, Development and Innovation Office -- NKFIH K124353 in the framework of ERA Chemistry. This material is based upon work supported by the National Science Foundation under Grant No. 1402897 (to D.G.). Present article was published in the frame of the project GINOP-2.3.2-15-2016-00053 (‘‘Development of engine fuels with high hydrogen content in their molecular structures (contribution to sustainable mobility)’’). Supported by the UNKP-17-4 (to MV) New National Excellence Program of the Ministry of Human Capacities. We appreciate the help of Adelina Voukadinova in thoroughly checking the database. We thank Zoltán Ható for his help with building the database.
-  G. Gouy. Sur la constitution de la charge electrique a la surface d’un electrolyte. J. Phys. (Paris), 9:457--468, 1910.
-  D. L. Chapman. A contribution to the theory of electrocapillarity. Phil. Mag., 25(148):475--481, 1913.
-  O. Stern. Zur theorie der elektrolytischen doppelschicht. Zeit. Elektrochemie, 30(21-22):508--516, 1924.
-  P. Simon and Y. Gogotsi. Materials for electrochemical capacitors. Nature Materials, 7(11):845--854, 2008.
-  S. Porada, R. Zhao, A. van der Wal, V. Presser, and P.M. Biesheuvel. Review on the science and technology of water desalination by capacitive deionization. Progress in Materials Science, 58(8):1388--1442, 2013.
-  H. U. Khan, M. E. Roberts, O. Johnson, R. Förch, W. Knoll, and Z. Bao. In situ, label-free DNA detection using organic transistor sensors. Advanced Materials, 22(40):4452--4456, 2010.
-  L. Friedrich and D. Gillespie. Improving charge-sensitive biomolecule sensors with the right choice of electrolyte. Sensors and Actuators B: Chemical, 230:281--288, 2016.
-  S. Haldrup, J. Catalano, M. Hinge, G. V. Jensen, J. S. Pedersen, and A. Bentien. Tailoring membrane nanostructure and charge density for high electrokinetic energy conversion efficiency. ACS Nano, 10(2):2415--2423, 2016.
-  F. H. J. van der Heyden, D. J. Bonthuis, D. Stein, C. Meyer, and C. Dekker. Electrokinetic energy conversion efficiency in nanofluidic channels. Nano Letters, 6(10):2232--2237, 2006.
-  F. H. J. van der Heyden, D. J. Bonthuis, D. Stein, C. Meyer, and C. Dekker. Power generation by pressure-driven transport of ions in nanofluidic channels. Nano Letters, 7(4):1022--1025, 2007.
-  D. Gillespie. High energy conversion efficiency in nanofluidic channels. Nano Letters, 12(3):1410--1416, 2012.
-  J.J. Bikerman. XXXIX. Structure and capacity of electrical double layer. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 33(220):384--397, 1942.
-  L. Blum and G. Stell. Solution of Ornstein-Zernike equation for wall-particle distribution function. J. Stat. Phys., 15:439--449, 1976.
-  L. Blum. Theory of electrified interfaces. J. Phys. Chem., 81(2):136--147, 1977.
-  D. Henderson, L. Blum, and W. R. Smith. Application of the hypernetted chain approximation to the electric double-layer at a charged planar interface. Chem. Phys. Lett., 63(2):381--383, 1979.
-  C.W. Outhwaite, L.B. Bhuiyan, and S. Levine. Theory of the electric double layer using a modified Poisson-Boltzman equation. J. Chem. Soc., Faraday Trans. 2, 76:1388--1408, 1980.
-  L.B. Bhuiyan, C.W. Outhwaite, and S. Levine. Numerical solution of a modified Poisson-Boltzmann equation for 1:2 and 2: 1 electrolytes in the diffuse layer. Mol. Phys., 42(6):1271--1290, 1981.
-  M. Lozada-Cassou, R. Saavedra-Barrera, and D. Henderson. The application of the hypernetted chain approximation to the electrical double-layer - comparison with Monte-Carlo results for symmetric salts. J. Chem. Phys., 77(10):5150--5156, 1982.
-  M. Lozada-Cassou and D. Henderson. Application of the hypernetted chain approximation to the electrical double-layer - comparison with Monte Carlo results for 2-1 and 1-2 salts. J. Phys. Chem., 87(15):2821--2824, 1983.
-  M. Plischke and D. Henderson. Pair correlation-functions and density profiles in the primitive model of the electric double-layer. J. Chem. Phys., 88(4):2712--2718, 1988.
-  Y. Rosenfeld. Free-energy model for the inhomogeneous hard-sphere fluid mixture and density-functional theory of freezing. Phys. Rev. Lett., 63(9):980--983, 1989.
-  E. Kierlik and M. L. Rosinberg. Free-energy density functional for the inhomogeneous hard-sphere fluid - application to interfacial adsorption. Phys. Rev. A, 42(6):3382--3387, 1990.
-  L. Mier y Teran, S. H. Suh, H. S. White, and H. T. Davis. A nonlocal free-energy density-functional approximation for the electrical double-layer. J. Chem. Phys., 92(8):5087--5098, 1990.
-  E. Kierlik and M. L. Rosinberg. Density-functional theory for inhomogeneous fluids - adsorption of binary-mixtures. Phys. Rev. A, 44(8):5025--5037, 1991.
-  Y. Rosenfeld. Free-energy model for inhomogeneous fluid mixtures - Yukawa-charged hard-spheres, general interactions, and plasmas. J. Chem. Phys., 98(10):8126--8148, 1993.
-  Y. Rosenfeld, M. Schmidt, H. Löwen, and P. Tarazona. Fundamental-measure free-energy density functional for hard spheres: Dimensional crossover and freezing. Phys. Rev. E, 55(4):4245--4263, 1997.
-  B. Groh, R. Evans, and S. Dietrich. Liquid-vapor interface of an ionic fluid. Phys. Rev. E, 57(6):6944--6954, 1998.
-  D. Boda, W. R. Fawcett, D. Henderson, and S. Sokołowski. Monte Carlo, density functional theory, and Poisson-Boltzmann theory study of the structure of an electrolyte near an electrode. J. Chem. Phys., 116(16):7170--7176, 2002.
-  D. Gillespie, W. Nonner, and R. S. Eisenberg. Coupling Poisson-Nernst-Planck and density functional theory to calculate ion flux. J. Phys.: Cond. Matt., 14(46):12129--12145, 2002.
-  D. Gillespie, W. Nonner, and R. S. Eisenberg. Density functional theory of charged, hard-sphere fluids. Phys. Rev. E, 68(3):031503, 2003.
-  D. di Caprio, J. Stafiej, and J. P. Badiali. Field theory for ionic systems. From fluctuations and structure at a hard wall to thermodynamics. Electrochim. Acta, 48(20-22):2967--2974, 2003.
-  D. di Caprio, J. Stafiej, and J. P. Badiali. Field theoretical approach to inhomogeneous ionic systems: thermodynamic consistency with the contact theorem, Gibbs adsorption and surface tension. Mol. Phys., 101(16):2545--2558, 2003.
-  O. Pizio, A. Patrykiejew, and S. Sokołowski. Phase behavior of ionic fluids in slitlike pores: A density functional approach for the restricted primitive model. J. Chem. Phys., 121(23):11957--11964, 2004.
-  D. Gillespie, M. Valiskó, and D. Boda. Density functional theory of the electrical double layer: the RFD functional. J. Phys.-Cond. Matt., 17(42):6609--6626, 2005.
-  J. Reszko-Zygmunt, S. Sokołowski, D. Henderson, and D. Boda. Temperature dependence of the double layer capacitance for the restricted primitive model of an electrolyte solution from a density functional approach. J. Chem. Phys., 122(8):084504, 2005.
-  M. Valiskó, D. Gillespie, and D. Boda. Selective adsorption of ions with different diameter and valence at highly-charged interfaces. J. Phys. Chem. C, 111(43):15575--15585, 2007.
-  T. Goel, C. N. Patra, S. K. Ghosh, and T. Mukherjee. A self-consistent density-functional approach to the structure of electric double layer: charge-asymmetric electrolytes. Mol. Phys., 107(1):19--25, 2009.
-  D. Henderson, S. Lamperski, Z. Jin, and J. Wu. Density functional study of the electric double layer formed by a high density electrolyte. J. Phys. Chem. B, 115(44):12911--12914, 2011.
-  M. Z. Bazant, B. D. Storey, and A. A. Kornyshev. Double layer in ionic liquids: Overscreening versus crowding. Phys. Rev. Lett., 106(4), 2011.
-  Z. Wang, L. Liu, and I. Neretnieks. The weighted correlation approach for density functional theory: a study on the structure of the electric double layer. J. Phys.: Cond. Matt., 23(17):175002, 2011.
-  O. Pizio, S. Sokołowski, and Z. Sokołowska. Electric double layer capacitance of restricted primitive model for an ionic fluid in slit-like nanopores: A density functional approach. J. Chem. Phys., 137(23):234705, 2012.
-  A. L. Frischknecht, D. O. Halligan, and M. L. Parks. Electrical double layers and differential capacitance in molten salts from density functional theory. J. Chem. Phys., 141(5):054708, 2014.
-  J. Jiang, D. Cao, D. Henderson, and J. Wu. Revisiting density functionals for the primitive model of electric double layers. J. Chem. Phys., 140(4):044714, 2014.
-  J. Jiang, D. Cao, D. Henderson, and J. Wu. A contact-corrected density functional theory for electrolytes at an interface. Phys. Chem. Chem. Phys., 16:3934--3938, 2014.
-  S. Zhou, S. Lamperski, and M. Zydorczak. Properties of a planar electric double layer under extreme conditions investigated by classical density functional theory and Monte Carlo simulations. J. Chem. Phys., 141(6):064701, 2014.
-  E.-Y. Kim, S.-C. Kim, Y.-S. Han, and B.-S. Seong. Structure of a planar electric double layer containing size-asymmetric ions: density functional approach. Mol. Phys., 113(8):871--879, 2015.
-  K. Ma, J. Forsman, and C. E. Woodward. Influence of ion pairing in ionic liquids on electrical double layer structures and surface force using classical density functional approach. J. Chem. Phys., 142(17):174704, 2015.
-  G. Yang and L. Liu. A systematic comparison of different approaches of density functional theory for the study of electrical double layers. J. Chem. Phys., 142(19):194110, 2015.
-  R. Roth and D. Gillespie. Shells of charge: a density functional theory for charged hard spheres. J. Phys.: Cond. Matt., 28(24):244006, 2016.
-  G. M. Torrie and J. P. Valleau. Electrical double-layers 4. Limitations of the Gouy-Chapman theory. J. Phys. Chem., 86(16):3251--3257, 1982.
-  D. Henderson and D. Boda. Insights from theory and simulation on the electrical double layer. Phys. Chem. Chem. Phys., 11(20):3822--3830, 2009.
-  G. M. Torrie and J. P. Valleau. Monte-Carlo study of an electrical double-layer. Chem. Phys. Lett., 65(2):343--346, 1979.
-  G. M. Torrie and J. P. Valleau. Electrical double-layers 1. Monte Carlo study of a uniformly charged surface. J. Chem. Phys., 73(11):5807--5816, 1980.
-  J. P. Valleau and G. M. Torrie. The electrical double-layer 3. Modified Gouy-Chapman theory with unequal ion sizes. J. Chem. Phys., 76(9):4623--4630, 1982.
-  J. P. Valleau and G. M. Torrie. Electrical double-layers 5. Asymmetric ion wall interactions. J. Chem. Phys., 81(12):6291--6295, 1984.
-  W. van Megen and I. Snook. The grand canonical ensemble Monte-Carlo method applied to the electrical double-layer. J. Chem. Phys., 73(9):4656--4662, 1980.
-  I. Snook and W. van Megen. Finite ion size effects in the electrical double-layer - a Monte-Carlo study. J. Chem. Phys., 75(8):4104--4106, 1981.
-  J. Loessberg-Zahl, K. G. H. Janssen, C. McCallum, D. Gillespie, and S. Pennathur. (Almost) stationary isotachophoretic concentration boundary in a nanofluidic channel using charge inversion. Anal. Chem., 88(12):6145--6150, 2016.
-  D. Gillespie. A review of steric interactions of ions: Why some theories succeed and others fail to account for ion size. Microfluidics and Nanofluidics, 18(5-6):717--738, 2014.
-  A. S. Khair and T. M. Squires. Ion steric effects on electrophoresis of a colloidal particle. J. Fluid Mech., 640:343--356, 2009.
-  E. Spohr. Some recent trends in computer simulations of aqueous double layers. Electrochim. Acta, 49(1):23--27, 2003.
-  E. Waisman and J. L. Lebowitz. Mean spherical model integral equation for charged hard spheres I. Method of solution. J. Chem. Phys., 56(6):3086--3093, 1972.
-  Y. He, D. Gillespie, D. Boda, I. Vlassiouk, R. S. Eisenberg, and Z. S. Siwy. Tuning transport properties of nanofluidic devices with local charge inversion. J. Am. Chem. Soc., 131(14):5194--5202, 2009.
-  D. Gillespie, A. S. Khair, J. P. Bardhan, and S. Pennathur. Efficiently accounting for ion correlations in electrokinetic nanofluidic devices using density functional theory. J. Coll. Interf. Sci., 359(2):520--529, 2011.
-  J. Hoffmann and D. Gillespie. Ion correlations in nanofluidic channels: Effects of ion size, valence, and concentration on voltage- and pressure-driven currents. Langmuir, 29(4):1303--1317, 2013.
-  D. Boda, K. Y. Chan, and D. Henderson. Monte Carlo simulation of an ion-dipole mixture as a model of an electrical double layer. J. Chem. Phys., 109(17):7362--7371, 1998.
-  A. Malasics, D. Gillespie, and D. Boda. Simulating prescribed particle densities in the grand canonical ensemble using iterative algorithms. J. Chem. Phys., 128(12):124102, 2008.
-  A. Malasics and D. Boda. An efficient iterative grand canonical Monte Carlo algorithm to determine individual ionic chemical potentials in electrolytes. J. Chem. Phys., 132(24):244103, 2010.
-  D. Boda and D. Gillespie. Calculating the electrostatic potential profiles of double layers from simulation ion density profiles. Hung. J. Ind. Chem., 41(2):125--132, 2013.