Microsoft Word - inner cover.docx HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(2) pp. 123–130 (2013) CALCULATING THE ELECTROSTATIC POTENTIAL PROFILES OF DOUBLE LAYERS FROM SIMULATION ION DENSITY PROFILES DEZSŐ BODA,B1 AND DIRK GILLESPIE2 1Department of Physical Chemistry, University of Pannonia, Egyetem u. 10., Veszprém, 8200, HUNGARY 2Department of Molecular Biophysics and Physiology, Rush University Medical Center, 1750 W. Harrison St., Chicago, USA BE-mail: boda@almos.vein.hu Computer simulations of the planar double layer geometry provide the charge profile with statistical noise. To compute the mean electrostatic potential profile from the charge profile, one must solve Poisson’s equation with appropriate boundary conditions (BC). In this work, we show that it is advantageous to use the Neumann or Dirichlet BCs at the boundaries of the simulation domain with an integrated version of Poisson’s equation. This minimises errors from the simulation’s noisy density profiles, in contrast to traditional convolution integrals that amplify the noise. The Neumann BC, where the electric field is prescribed, can be used in both the constant surface charge and constant electrode voltage ensembles. In the constant voltage ensemble, where the potential difference between the confining electrodes is prescribed, one can also use the Dirichlet BC, where the potentials at the boundaries are set. We show that the new methods provide converged results for the potential profile faster than the convolution integral does. Keywords: electrical double layer, Poisson equation, boundary conditions, computer simulation Introduction The electrical double layer (DL) is formed by a charged surface and a phase (usually, a liquid) containing mobile charge carriers near the surface. Depending on the ma- terial carrying the mobile charges, DLs appear in elec- trolytes, molten salts, ionic liquids, plasmas, and even fast ion conductors (solid electrolytes). DLs in solutions of dissolved ions are particularly important in electrochem- istry, biology, and colloid chemistry. DLs near electrodes differ from DLs near charged objects carrying a fixed surface charge (such as colloids, macromolecules, and porous bodies), because the surface charge on the elec- trode can be controlled by an external voltage. Theoretical studies of DLs began with the Poisson- Boltzmann (PB) theory known as the Gouy-Chapman (GC) theory [1, 2] in electrochemistry, the Debye-Hückel (DH) theory [3] in solution chemistry, and the Derjaguin- Landau-Verwey-Overbeek (DLVO) theory [4, 5] in col- loid chemistry. The PB theory is still very popular in ap- plications because of its simplicity in spite of the fact that it neglects ionic correlations and effects due to the fi- nite size of ions (e.g., excluded volume). More powerful statistical mechanical theories that are able to take these correlations into account have been developed [6–16]. Computer simulations are a versatile method of study- ing DLs in various geometries for various models of the constituents (ions, water, and electrode) [17–35]. This pa- per focuses on computing the electrical field and potential from a charge distribution obtained from computer simu- lations. Computer simulations of systems containing charged particles must be performed in accordance with the laws of both electrostatics and statistical mechanics [36, 37]. This means that the electrical potential must be com- puted accurately for every configuration sampled in a Monte Carlo (MC) simulation. Alternatively, the electri- cal field must be computed accurately in every time step of a molecular dynamics (MD) simulation. Every con- figuration of the system (positions of ions) corresponds to a microscopic state. In this case, we consider the sys- tem at the microscopic level. A simulation that samples the possible microscopic states must be performed prop- erly according to the probability distribution of a given statistical mechanical ensemble. One typical simulation method to handle various ensembles is MC. The density distribution of ionic species i, h⇢ i (r)i, is obtained as an ensemble average from the simulation. From these, the average charge distribution can be obtained as q(r) = hQ(r)i = X i z i e h⇢ i (r)i , (1) where Q(r) is the charge distribution in a microscopic state, e is the electronic unit charge and z i is the valence of ionic species i. The average (mean) electrostatic potential can be ob- tained “on the fly” by computing the potential in the sim- ulation cell for every configuration (denoted by (r)) and then taking the ensemble average, (r) = h (r)i. In this work, we use upper-case symbols for the micro- scopic quantities ( (r) and Q(r)), while we use lower- case symbols for their macroscopic counterparts, namely, 124 their ensemble averages ( (r) and q(r)). The electrical potential in a configuration of the MC simulation, (r), can be computed analytically from Coulomb’s law or nu- merically using a Poisson solver. In both cases, we apply electrostatic BCs at the microscopic level. The reverse order, when we compute the ensemble av- erage of the charge distribution (Eq.(1)) and then solve Poisson’s equation for the mean potential r2 (r) = � 1 ✏ 0 q(r) (2) is more usual (✏ 0 is the permittivity of vacuum). This equation applies to either explicit or implicit solvent models. In the explicit solvent framework, q(r) also con- tains the charge distribution of the water molecules (in addition to their ionic charges). In the implicit solvent framework, where water is represented by a dielectric background characterised by a dielectric constant, ✏, q(r) also contains polarisation charges induced in the dielec- tric (in the simplest case, q(r) is the ionic charges divided by ✏). Because we are solving a differential equation, we must use appropriate boundary conditions (BCs) when solving Eq.(2). Because we want an ensemble averaged result, we apply the BCs at the macroscopic level. This paper describes how to numerically integrate the Poisson equation with appropriate BCs to efficiently compute the mean electrostatic potential profile from the mean charge profile obtained from simulations. The sta- tistical ensemble applied in the simulation determines which BC to use. In this paper, we consider the two ba- sic ensembles, one where the electrode charges are fixed (constant charge ensemble) and one where the difference between the electrode potentials is fixed (constant voltage ensemble). In the constant charge ensemble, the traditional method of computing the potential profile is via a convo- lution integral. Here, we show that this method is error- prone and numerically inefficient and instead, it is more advantageous to use a different integration scheme with Neumann BCs, where the normal electric fields at the boundaries of the system are fixed [38]. In the constant voltage ensemble developed by KIY- OHARA and ASAKA [39], the electrode potentials are known in advance. Therefore, we can also use Dirichlet BCs, where the potentials at the boundaries of the system are prescribed. In the following, we describe our model in detail. Then, we consider all three issues discussed above (BCs at the microscopic level, statistical mechanical ensem- bles, and BCs at the macroscopic level) and present var- ious possibilities for the macroscopic BC depending on the statistical mechanical ensemble used. We present re- sults of model calculations to show the self-consistency of these calculations. Model and Boundary Condition at the Microscopic Level At the microscopic level, we have a system that con- tains localised discrete and/or continuous distributions of charges Q(r). In practice, there are two traditional schools to compute the electrostatic energy (in MC) or forces (in MD) for a configuration sampled by a com- puter simulation. In one school, the potential is computed on a grid from Poisson’s equation (Eq.(2)) using a par- tial differential equation (PDE) solver with appropriate boundary conditions. This method is generally used in MD simulations of explicit solvent systems. In the other school, Coulomb’s law is used to calculate the potential: (r) = 1 4⇡✏ 0 ✏ Z Q(r0)G(r, r0)dr0, (3) where G(r, r0) is the appropriate Green’s function. This method is generally used for simulations with implicit solvents. The planar DL geometry means that we have a rect- angular simulation cell of length L 2 �L 1 and with a base H ⇥ H confined by two planar electrodes at the two ends (x = L 1 and x = L 2 ) carrying surface charges � 1 and � 2 . Periodic boundary conditions (PBC) are used in the y, and z dimensions, which means that the Green’s func- tion is G(r, r0) = 1X j=�1 1X k=�1 1 |r � r0 + jHn y + kHn z | (4) where the sum over j and k represents the interaction with the periodic image charges in the replicas of the cen- tral simulation cell in the y and z dimensions set by the unit vectors n y and n z . The term j = k = 0 corresponds to the interaction with the charge in the central simula- tion cell computed explicitly. In this paper, the interac- tion with periodic replicas is taken into account using the charged sheet method by smearing these charges into a sheet carrying q/H2 surface charges with a H2 square hole in the middle. The interaction with the holed sheet can be integrated. For further details, see Ref. [18, 29]. The statistical mechanical ensemble determines which thermodynamic variables are fixed in the simula- tion. The attempts in the MC simulations are designed to ensure sampling according to the probability distribution of the given ensemble. Moreover, the simulation must be self-consistent in the sense that a prescribed thermody- namic variable must agree with its value computed as an output of the simulation. For example, in an NpT simu- lation (where the pressure is fixed), the pressure can also be computed from the virial sum as an ensemble aver- age in both NV T and NpT simulations. This value must be equal to the one prescribed in the NpT ensemble. In practice, however, the accuracy of the pressure computed from the virial sum depends on the size of the system, and one obtains the same relationship between pressure 125 and density from NV T and NpT simulations only in the limit of very large simulation cells. From the point of view of the DL problem, the ques- tion is whether we perform the simulation in the con- stant charge or constant voltage ensemble. In the con- stant charge ensemble the surface charges on the con- fining walls of the simulation cell are fixed. This is the traditional simulation setup used in DL simulations since TORRIE and VALLEAU [17, 18]. The output of the simu- lation is the density profiles of the various ionic species h⇢ i (x)i (the PBC applied in the y and z dimensions en- sures that the profiles depend on x only), from which the x-dependent charge profile, q(x), is obtained (see Eq.(1)). The corresponding form of Poisson’s equation is d2 (x) dx2 = � 1 ✏ 0 ✏ q(x), (5) where ✏ is now included so q(x) constrains only ionic and electrode charges. Recently, KIYOHARA and ASAKA [39] introduced the constant voltage ensemble, where the potential differ- ence between the two confining walls (the electrodes) is prescribed. A special MC step, where a small amount of charge, ��, is transferred from one electrode to the other, has been introduced. The charge exchange is accepted on the basis of the Boltzmann factor min  1, exp ✓ � �U kT + H2� ⇤�� kT ◆� (6) where �U is the energy change associated with the charge movement and � ⇤ = (L 2 ) � (L 2 ) is the prescribed voltage. Here, the surface charges on the elec- trodes fluctuate, while the potential difference between the electrodes (� ⇤) is an independent, prescribed vari- able of the ensemble. The mean potential profile can be computed from Eq.(5). The computed voltage, namely, the potential difference between the electrodes, � = (L 2 ) � (L 2 ), must be equal to the prescribed voltage, � ⇤. In the following, we describe various ways to apply BCs at the macroscopic level using a case study of a 1:1 electrolyte, where the ions are modeled as charged hard spheres with diameters d + = d� = 3 Å and the di- electric constant of water is 78.46. The concentration is 1 M and the temperature is 298.2 K. We show results for the special case of � 1 = �� 2 = � = 0.1 Cm�2 and L 1 = �L 2 = L, but the equations are presented for the general case. We used 200/200 ions in the MC simula- tions performed in the canonical (NV T ) ensemble. The dimensions of the cell are L = 50 Å and H = 57.2 Å. Boundary Conditions at the Macroscopic Level Boundary Conditions Set in the Bath: The Convolution Integral Let us distinguish between the charge of the ions obtained from the simulation, q ion (x), and the surface charges, P k � k �(x � x k ), at x 1 = L 1 and x 2 = L 2 (in this work, we assume only two charged surfaces at L 1 and L 2 , but there can be more as in the simulations of KIYOHARA et al. [40–50] for porous electrodes). The total charge is the sum of these: q(x) = q ion (x) + X k � k �(x � x k ). (7) The ionic charge profile is obtained as an ensemble aver- age. The electrode charges are prescribed in the constant charge ensemble, while they are obtained as ensemble av- erages in the constant voltage ensemble. At the macroscopic level, the issue of electrostatic self-consistency appears when we ask the question: what kind of BC should be applied when we solve Poisson’s equation (Eq.(5)). The traditional answer to this question is that the BC is set in the bulk electrolyte, where the av- erage electrical field and electrical potential are zero. The corresponding solution of Poisson’s equation can then be obtained in the form of a convolution integral (x) = � 1 ✏ 0 ✏ Z 1 x (x0 � x)q ion (x0)dx0. (8) This solution was probably inspired by theories that usually consider an isolated DL where the BCs are set in infinity. Because theories (unlike simulations) provide smooth charge profiles (with the property lim x!1 q(x) = 0) without any noise, this integral works well for theories. To the best of our knowledge, the ma- jority of researchers (among others, the authors of this paper) have used this equation in the past [6–35]. In this work, we show that this equation, from a numerical point of view, is a poor choice to compute the potential in sim- ulation studies. The reason is that the upper integration limit is not well defined and that q(x) is subject to a large statistical noise. In practice, the upper integration limit is set some- where in the middle of the cell where a bulk electrolyte is. Here we will use exactly the middle of the cell (x = 0) as the upper integration limit. Because of x�x0 in the in- tegrand, this integral is very sensitive to the noise in q(x), because the noise is amplified as one moves further away from the electrode. The results of a very short (200 MC cycles; 1000 at- tempts to move ions were made in an MC cycle) simu- lation are shown in Fig.1. The density and charge pro- files are very noisy. When we compute the potential from Eq.(8), the result is subject to a large error and is far from what we are supposed to get (Fig.2, top panel). The slope is not necessarily zero in the bath (which means that the electrical field is not zero). If the simulation was run for longer, we would get a completely different re- sult. In general, long simulations are needed to produce a smooth charge profile and a well established potential profile. The problem is even more serious when we try to reproduce small effects, such as the value of the electrode potential at zero electrode charge (PZC) for asymmetric electrolytes. The value of the PZC potential is very small 126 0 1 2 3 4 5 D en si ty p ro fi le s / M -60 -40 -20 0 20 40 60 x / Å -4 -2 0 2 4 C ha rg e pr of il e / M Figure 1: Density profiles (top panel) and a charge profile (bottom panel) obtained from a short (200 MC cycles) constant charge simulation using � 1 = �� 2 = �0.1 Cm�2. The charge profile is obtained from P i z i ⇢ i (x), so its unit is M (mol dm�3) and the effect of the noise is dramatic. Extremely long simulations are needed to obtain convergent results for the potential [51]. To illustrate this weak convergence, we have plotted the left electrode potential ( (L 1 )� (0), top panel) and the voltage ( (L 2 ) � (L 1 ), bottom panel) as computed from Eq.(8) as functions of the performed MC cycles in Fig.3 (red dashed curves). These potential values fluctu- ate strongly. The calculated values depend not only on the simulation time, but on the upper integration limit. If we shift that point a bit, we get a different result (data not shown). Neumann Boundary Conditions and the Constant Charge Ensemble Here we propose, instead, to use Neumann or Dirichlet BCs at the boundaries of the simulation cell (at L 1 and L 2 , or, equivalently, at �1 and 1). In the case of the Neumann BC, the normal electrical field is prescribed, while in the case of the Dirichlet BC, the electrical poten- tial is prescribed at the confining walls. Our simulation setup ensures that the simulation cell is always charge neutral. Then, Gauss law states that the average electrical field is zero outside the cell for the regions z < L 1 and z > L 2 . This information makes the Neumann BC ap- plicable both in the constant charge and constant voltage 0 1 2 3 4 ψ (x ) Neumann convolution (left) convolution (right) -60 -40 -20 0 20 40 60 x / Å 0 0.1 0.2 0.3 0.4 0.5 0.6 q t ot (x ) Figure 2: Potential profiles (top panel) as computed from the Neumann BC (Eq.(11)) and the convolution integral (Eq.(8)). The upper limit of the convolution integral is x = 0 and it is computed for the left- and right-hand sides separately. The bottom panel shows the integral of the charge profile (q tot (x), see Eq.(7)). It is closely related to the electric field through Eq.(13). The profiles have been obtained from the curves of Fig.1. The electrostatic potential is shown in units of kT/e throughout this paper ensembles. In the constant charge ensemble, the potential differ- ence between the electrodes is an output of the calcula- tion. Therefore, we cannot use it as a BC, so we cannot apply the Dirichlet BC in this case. In the constant volt- age ensemble, on the other hand, the potential difference is known in advance (see the next section). In the constant voltage ensemble, therefore, both Neumann and Dirichlet BCs can be used and they should give the same answer (apart from errors related to the size of the system, see later). The Neumann BC for Eq.(8) is that the electrical field is zero outside the system: E(x ! �1) = � d (x) dx ���� x!�1 = 0 (9) and the same for x ! 1. By integrating Poisson’s equa- tion once, we obtain d (x) dx = � 1 ✏ 0 ✏ Z x �1 q(x0)dx0 + C 1 (10) where C 1 is an integration constant. Taking Eq.(10) at 127 -4 -3.5 -3 -2.5 -2 -1.5 -1 ψ (L 1) -ψ (0 ) Neumann convolution 0 200 400 600 800 1000 MC cycle 2 3 4 5 6 ψ (L 2) -ψ (L 1) Figure 3: Convergence of the left electrode potential (top panel) and the voltage (bottom panel) as computed from the Neumann BC (Eq.(11)) and the convolution integral (Eq.(8)) in the constant charge ensemble. any location x < L 1 (where q(x) = 0) and using the BC (Eq.(9)), we get C 1 = 0 for the integration constant. If we use Eq.(7) for q(x), we obtain d (x) dx = � 1 ✏ 0 ✏ Z x L1 q ion (x0)dx0 � 1 ✏ 0 ✏ � 1 (11) for L 1 < x < L 2 . The change in the lower integration limit was possible because ions exist only between the two electrodes (q ion 6= 0 only for L 1 < x < L 2 ). Intro- ducing q tot (x) = Z x L1 q ion (x0)dx0 (12) for the integral of the ionic charge profile (the total charge density per area in the [L 1 , x] interval), the electrical field can be given as E(x) = 1 ✏ 0 ✏ q tot (x) + 1 ✏ 0 ✏ � 1 (13) for L 1 < x < L 2 . The q tot (x) profile is shown in the bottom panel of Fig.2. For the special case of � 1 = �� 2 , the total ionic charge is zero in the system, so q tot (L 2 ) = 0. The nearly constant but noisy profile in the middle of the cell represents the bulk region. Integrating once more, we obtain (x) = � 1 ✏ 0 ✏ Z x L1 q tot (x0)dx0 � 1 ✏ 0 ✏ � 1 (x � L 1 ) + C 2 , (14) -3.5 -3 -2.5 -2 ψ (L 1) -ψ (0 ) Dirichlet Neumann convolution 0 200 400 600 800 1000 MC cycle 4 5 6 7 ψ (L 2) -ψ (L 1) Figure 4: Convergence of the left electrode potential (top panel) and the voltage (bottom panel) as computed from the Neumann BC (Eq.(11)), the Dirichlet BC (Eq.(17)), and the convolution integral (Eq.(8)) in the constant voltage ensemble where C 2 is another integration constant. Its value is in- consequential because we can set the zero level of the potential arbitrarily. Fig.2 shows the (x) profile where C 2 = 0. It is usual, however, to set the ground in the bulk, so C 2 = � (0), where (0) is the average of the potential profile over the bulk. KIYOHARA and ASAKA [39, 52] used an equation of a convolution form that can be shown to be equivalent to Eq.(14), but they seemed to leave out the linear term containing � 1 . Eq.(14) is different from the convolution integral (Eq.(8)) because it does not suffer from the uncertainty in the upper integration limit and from errors originating from the noise of the q ion (x) profile. The integration is performed for a well defined finite domain from the left electrode to x. Fig.3 shows that this equation provides a much better convergence as a function of simulation time for the same simulation (same q ion (x)). This procedure was used in our papers for inhomoge- neous electrolyte systems to compute the mean electro- static potential [53–56]. It was especially useful for elec- trolytes adsorbed in narrow slits [55, 56]. In this case, the DLs formed at the walls of the slit overlap so a bulk elec- trolyte does not form in the middle of the slit and the slit is not charge neutral. The Neumann BC is then the nat- ural BC so the electric field is zero behind the walls. In another example, the Neumann BC is used to compute the potential for a DL model, where the electrode, the 128 inner layer, and the electrolyte have different dielectric constants [54]. In this case, the polarisation charge in- duced at the dielectric boundaries must also be included in Eq.(14). Dirichlet Boundary Conditions and the Constant Voltage Ensemble An alternative method is the constant voltage ensem- ble of KIYOHARA and ASAKA [39]. Here, the electrode charges fluctuate and the potential between the electrodes is prescribed. Therefore, we can also apply Dirichlet BCs, where (L 1 ) = 0 (15) and (L 2 ) = � ⇤. (16) The general solution for the potential profile in the L 1 < x < L 2 range is (x) = � 1 ✏ 0 ✏ Z x L1 q tot (x0)dx0 +C 1 (x�L 1 )+C 2 , (17) where we integrate from the right-hand side of the elec- trode at x = L 1 , so the surface charge � 1 is now excluded from the integration. The BC at x = L 1 (Eq.(15)) pro- vides the integration constant C 2 = 0, while the BC at x = L 2 (Eq.(16)) provides the integration constant C 1 = 1 L 2 � L 1 " � ⇤ + 1 ✏ 0 ✏ Z L2 L1 q tot (x)dx # . (18) Of course, we can calculate the potential profile using the Neumann BCs too. In that case, we must use h� 1 i in Eq.(14) instead of � 1 because the electrode charge is now fluctuating so its value is not known in advance. There- fore, its ensemble average should be used in Eqs.(7-14). In the constant voltage ensemble we need the value of the voltage that corresponds to � 1 = �0.1 Cm�2 as used in the previous constant charge simulation. This value we estimated with a very long (50,000 MC cy- cles) constant charge simulation and was obtained as � ⇤ = 4.398 kT/e. This value was used in the constant voltage simulation as an input parameter. Fig.4 is the analogous version of Fig.3. The black solid curves are the results obtained from the Dirichlet BCs (Eqs.(17) and (18)). The red short-dashed curves show the results of the Neumann BC using h� 1 i in Eq.(14). The convolution integral results (blue long- dashed curves) are inaccurate and poorly converged for such a short simulation. In the constant voltage ensemble the electrode charge is a fluctuating quantity. Fig.5 shows its convergence. Its limiting value is not equal to that used in the constant charge simulation that provided the input voltage value � ⇤ = 4.398 kT/e. The deviation is due to finite size of the system. Using a larger simulation box (larger H), a smaller deviation is observed (data not shown). This de- viation is also observed in Fig.4. The limiting value of the 0 200 400 600 800 1000 MC cycle -0.1 -0.095 -0.09 E le ct ro de c ha rg e / C m -2 Figure 5: Convergence of the electrode charge in the constant voltage ensemble potential difference � = (L 2 )� (L 1 ) obtained from the calculation using the Neumann BC (bottom panel) is different from the prescribed value (the value used in the Dirichlet BC, see also bottom panel). This behaviour is analogous to that discussed earlier regarding the example of the pressure computed in the NV T and NpT ensembles. The constant charge ensem- ble corresponds to the NV T , while the constant volt- age ensemble corresponds to the NpT ensemble. In the constant charge ensemble, the voltage is computed us- ing the Neumann BC, just as the pressure is computed in the NV T ensemble from the virial sum. In the constant voltage ensemble, the voltage is prescribed, just as the pressure is prescribed in the NpT ensemble. The voltage can also be computed from the Neumann BC, just as the pressure can also be computed from the virial sum in the NpT ensemble. Conclusion We propose that the Neumann or Dirichlet BCs should be used in computing the mean electrostatic potential for the planar DL geometry studied by computer simulations. The commonly used convolution integral of Eq.(8) re- quires a vaguely defined upper integration limit and also suffers from numerical problems because it magnifies the effect of the noise in the charge profile that is always present in computer simulations. The problem of noise in the density profiles is unique to simulations and there- fore more care must be taken in computing the electro- static potential. On the other hand, theories [6–16], which produce smooth, noise-free density profiles, can use the convolution integral. We have shown here that the numerical method we proposed is much more efficient than the convolution in- tegral, because we use unambiguous parameters in the BC, specifically, the electrode charge in the case of the Neumann BC and the voltage in the case of the Dirichlet BC. Also, the simulation cell is necessarily finite, there- fore, the boundaries of the system are always well de- 129 fined. Overall, our method leads to converged results with very short simulations. Acknowledgement The present publication was realised with the support of the projects TÁMOP-4.2.2/A-11/1/KONV-2012-0071 and TÁMOP-4.1.1/C-12/1/KONV-2012-0017. REFERENCES [1] GOUY G.: Sur la constitution de la charge elec- trique a la surface d’un electrolyte (English Title Please), J. Phys. (Paris), 1910, 9, 457 (in French) [2] CHAPMAN D.L.: A Contribution to the Theory of Electrocapillarity, Phil. Mag., 1913, 25, 475 [3] DEBYE P., HÜCKEL E.: The theory of electrolytes. I. Lowering of freezing point and related phenom- ena, Physik. Z., 1923, 24, 185–206 [4] DERJAGUIN B., LANDAU L.: Theory of the stabil- ity of strongly charged lyophobic sols and of the ad- hesion of strongly charged particles in the solution of electrolytes, Acta Physiochem USSR, 1941, 14, 633–662 [5] VERWEY E.J.W., OVERBEEK J.T.G.: Theory of the Stability of Lyophobic Colloids Elsevier, Am- sterdam, 1948 [6] BLUM L.: Theory of electrified interfaces, J. Phys. Chem., 1977, 81(2), 136–147 [7] HENDERSON D., BLUM L., SMITH W.R.: Applica- tion of the hypernetted chain approximation to the electric double-layer at a charged planar interface, Chem. Phys. Lett., 1979, 63(2), 381–383 [8] BLUM L., HENDERSON D.: Mixtures of hard ions and dipoles against a charged wall: The Ornstein– Zernike equation, some exact results, and the mean spherical approximation, J. Chem. Phys., 1981, 74(3), 1902–1910 [9] OUTHWAITE C.W., BHUIYAN L.B.: An improved modified Poisson-Boltzmann equation in electric- double-layer theory, J. Chem. Soc. Faraday. Trans. II., 1983, 79, 707–718 [10] LOZADA-CASSOU M., HENDERSON D.: Applica- tion 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., 1983, 87(15), 2821–2824 [11] MIER-Y TERAN L., SUH S.H., WHITE H.S., DAVIS H.T.: A nonlocal free-energy density- functional approximation for the electrical double- layer, J. Chem. Phys., 1990, 92(8), 5087–5098 [12] KIERLIK E., ROSINBERG M.L.: Free-energy den- sity functional for the inhomogeneous hard-sphere fluid - application to interfacial adsorption, Phys. Rev. A, 1990, 42(6), 3382–3387 [13] ROSENFELD Y.: Free-energy model for inhomo- geneous fluid mixtures - Yukawa-charged hard- spheres, general interactions, and plasmas, J. Chem. Phys., 1993, 98(10), 8126–8148 [14] GILLESPIE D., NONNER W., EISENBERG R.S.: Density functional theory of charged, hard-sphere fluids, Phys. Rev. E, 2003, 68(3), 031503 [15] DI CAPRIO D., STAFIEJ J., BADIALI J.P.: Field theoretical approach to inhomogeneous ionic sys- tems: thermodynamic consistency with the con- tact theorem, Gibbs adsorption and surface tension, Mol. Phys., 2003, 101(16), 2545–2558 [16] PIZIO O., PATRYKIEJEW A., SOKOLOWSKI S.: Phase behavior of ionic fluids in slitlike pores: A density functional approach for the restricted prim- itive model, J. Chem. Phys., 2004, 121(23), 11957– 11964 [17] TORRIE G.M., VALLEAU J.P.: Monte Carlo study of an electrical double-layer, Chem. Phys. Lett., 1979, 65(2), 343–346 [18] TORRIE G.M., VALLEAU J.P.: Electrical double- layers 1. Monte Carlo study of a uniformly charged surface, J. Chem. Phys., 1980, 73(11), 5807–5816 [19] TORRIE G.M., VALLEAU J.P., PATEY, G.N.: Elec- trical double-layers 2. Monte Carlo and HNC stud- ies of image effects, J. Chem. Phys., 1982, 76(9), 4615–4622 [20] VALLEAU J.P., TORRIE G.M.: The electrical double-layer 3. Modified Gouy-Chapman theory with unequal ion sizes, J. Chem. Phys., 1982, 76(9), 4623–4630 [21] TORRIE G.M., VALLEAU J.P.: Electrical double- layers 4. Limitations of the Gouy-Chapman theory, J. Phys. Chem., 1982, 86(16), 3251–3257 [22] VALLEAU J.P., TORRIE G.M.: Electrical double- layers 5. Asymmetric ion wall interactions, J. Chem. Phys., 1984, 81(12), 6291–6295 [23] TORRIE G.M., VALLEAU J.P., OUTHWAITE C.W.: Electrical double-layers 6. Image effects for diva- lent ions, J. Chem. Phys., 1984, 81(12), 6296–6300 [24] VAN MEGEN W., SNOOK I.: The grand canoni- cal ensemble Monte Carlo method applied to the electrical double-layer, J. Chem. Phys., 1980, 73(9), 4656–4662 [25] SNOOK I., VAN MEGEN W.: Finite ion size effects in the electrical double-layer - a Monte Carlo study, J. Chem. Phys., 1981, 75(8), 4104–4106 [26] BRATKO D., JÖNSSON B., WENNERSTRÖM H.: Electrical double layer interactions with image charges, Chem. Phys. Lett., 1986, 128, 449–454 [27] TANG Z.X., SCRIVEN L.E., DAVIS H.T.: A 3- component model of the electrical double-layer, J. Chem. Phys., 1992, 97(1), 494–503 [28] PHILPOTT M.R., GLOSLI J.N.: Molecular dynam- ics simulation of interfacial electrochemical pro- cesses: Electric double layer screening, Solid-liquid Electrochemical Interfaces, 1997, 656, 13–30 130 [29] BODA D., CHAN K.Y., HENDERSON D.: Monte Carlo simulation of an ion-dipole mixture as a model of an electrical double layer, J. Chem. Phys., 1998, 109(17), 7362–7371 [30] SPOHR E.: Computer simulation of the structure of the electrochemical double layer, J. Electroanalyti- cal Chem., 1998, 450(2), 327–334 [31] BODA D., HENDERSON D., CHAN K.Y.: Monte Carlo study of the capacitance of the double layer in a model molten salt, J. Chem. Phys., 1999, 110(11), 5346–5350 [32] CROZIER P.S., ROWLEY R.L., HENDERSON D.: Molecular-dynamics simulations of ion size effects on the fluid structure of aqueous electrolyte systems between charged model electrodes, J. Chem. Phys., 2001, 114(17), 7513–7517 [33] GUYMON C.G., HUNSAKER M.L., HARB J.N., HENDERSON D., ROWLEY R.L.: Effects of solvent model flexibility on aqueous electrolyte behavior between electrodes, J. Chem. Phys., 2003, 118(22), 10195–10202 [34] SPOHR E.: Some recent trends in computer sim- ulations of aqueous double layers, Electrochimica Acta, 2003, 49(1), 23–27 [35] HENDERSON D., BODA D.: Insights from theory and simulation on the electrical double layer, Phys. Chem. Chem. Phys., 2009, 11(20), 3822–3830 [36] ALLEN M.P., TILDESLEY D.J.: Computer Simula- tion of Liquids Oxford, New York, 1987 [37] FRENKEL D., SMIT B.: Understanding molecular simulations Academic Press, San Diego, 1996 [38] JACKSON J.D.: Classical Electrodynamics Wiley, New York, 3rd edn., 1999 [39] KIYOHARA K., ASAKA K.: Monte Carlo simu- lation of electrolytes in the constant voltage en- semble, J. Chem. Phys., 2007, 126(21), 214704 (pages 14) [40] KIYOHARA K., ASAKA K.: Monte Carlo simula- tion of electrolytes in the constant voltage ensem- ble, J. Chem. Phys., 2007, 126(21), 214704 [41] KIYOHARA, K., ASAKA, K.: Monte Carlo simu- lation of porous electrodes in the constant voltage ensemble, J. Phys. Chem. C, 2007, 111(43), 5903– 15909 [42] KIYOHARA K., SUGINO T., TAKEUCHI I., MUKAI K., ASAKA K.: Expansion and contraction of poly- mer electrodes under an applied voltage, J. Appl. Phys., 2009, 105(6), 063506 [43] KIYOHARA K., SUGINO T., TAKEUCHI I., MUKAI K., ASAKA K.: Expansion and contraction of poly- mer electrodes under an applied voltage, J. Appl. Phys., 2009, 105(11), 119902 [44] KIYOHARA K., SUGINO T., ASAKA K.: Elec- trolytes in porous electrodes: Effects of the pore size and the dielectric constant of the medium, J. Chem. Phys., 2010, 132(14), 144705 [45] KIYOHARA K., SUGINO T., ASAKA K.: Phase transition in porous electrodes, J. Chem. Phys., 2011, 134(15), 154710 [46] KIYOHARA K., SUGINO T., ASAKA K.: Molecu- lar mechanism of ionic electroactive polymer actu- ators, Smart Mat. & Struc., 2011, 20(12), 124009 [47] KIYOHARA K., SHIOYAMA H., SUGINO T., ASAKA K.: Phase transition in porous electrodes. II. Effect of asymmetry on ion size, J. Chem. Phys., 2012, 136(9), 094701 [48] KIYOHARA K., ASAKA K.: Voltage induced pres- sure in porous electrodes, Mol. Phys., 2013, 111(2), 295–306 [49] KIYOHARA K., SHIOYAMA H., SUGINO T., ASAKA K., SONEDA Y., IMOTO K., KODAMA M.: Phase transition in porous electrodes. III. In the case of a two component electrolyte, J. Chem. Phys., 2013, 138(23), 234704 [50] KIYOHARA K., SHIOYAMA H., SUGINO T., ASAKA K., SONEDA Y., IMOTO K., KODAMA M.: Phase transition in porous electrodes. III. In the case of a two component electrolyte (vol 138, 234704, 2013), J. Chem. Phys., 2013, 139(6), 069902 [51] VALISKÓ M., HENDERSON D., BODA D.: Com- petition between the effects of asymmetries in ion diameters and charges in an electrical double layer studied by Monte Carlo simulations and theories, J. Phys. Chem. B, 2004, 108(42), 16548–16555 [52] KIYOHARA K., ASAKA K.: Monte Carlo Simula- tion of Porous Electrodes in the Constant Voltage Ensemble, J. Phys. Chem. C, 2007, 111(43), 15903– 15909 [53] NAGY T., VALISKÓ M., HENDERSON D., BODA D.: The Behavior of 2:1 and 3:1 Electrolytes at Polarizable Interfaces, J. Chem. Eng. Data, 2011, 56(4), 1316–1322 [54] NAGY T., HENDERSON D., BODA D.: Simulation of an electrical double layer model with a low di- electric layer between the electrode and the elec- trolyte, J. Phys. Chem. B, 2011, 115(39), 11409– 11419 [55] KOVÁCS R., VALISKÓ M., BODA D.: Monte Carlo simulation of the electrical properties of elec- trolytes adsorbed in charged slit-systems, Cond. Matt. Phys., 2012, 15(2), 23803 [56] VALISKÓ M., HENDERSON D., BODA D.: Se- lective adsorption of ions in charged slit systems, Cond. Matt. Phys., 2013, 16(4), 43601