Grand canonical simulations of ions between charged conducting surfaces using exact 3D Ewald summations

Abstract

We present a useful methodology to simulate ionic fluids confined by two charged and perfectly conducting surfaces. Electrostatic interactions are treated using a modified 3D Ewald sum, which accounts for all image charges across the conductors, as well as the 2D periodicity, parallel to the surfaces. The energy expression is exact, and the method is trivial to implement in existing Ewald codes. We furthermore invoke a grand canonical scheme that utilizes a bias potential, that regulates the surface charge density. The applied bias potential also enables us to calculate individual chemical potentials of the ions. Finally, we argue that our approach leads to a pedagogically appealing description of the Donnan potential, and what it measures in these systems.

---

1 Introduction

Mechanisms underlying different electrochemical properties of ionic systems confined by polarizable surfaces has attracted much attention recently due to their application in electric double layers capacitors.^1–3 Simulating such systems provides a number of challenges, associated mainly with treating the long-range electrostatic interactions that are used to obtain electrochemical properties such as surface potentials and capacitance. Many different methods have been developed to manage electrostatic interactions in ionic fluids confined between two flat conducting electrodes. Examples include solutions to the Poisson equation using expansion methods,⁴ the ICMMM2D method,⁵ the ELCIC method,⁶ mapping schemes,⁷ direct summations of the image charges^8,9 and Gaussian representation of the surface charges.¹⁰ However, generally such methods can be difficult to implement with slow convergence. Instead, in this work we implement a different approach which utilizes a double cell that includes the ions and their nearest images to simulate the many-reflected full system via Ewald summation. The Ewald summation has a long history in the simulation of coulombic fluids and has been developed over many decades to become an efficient and valuable tool for researchers. As such, optimized implementations, such as particle mesh Ewald¹¹ and smooth particle mesh Ewald¹² can be put to good use to efficiently treat the problem of charged fluids confined by conducting surfaces.

One of the difficulties associated with calculating properties such as the capacitance in these systems is that knowledge of the electrode potential relative to a fixed reference is required. While the latter is often chosen to be the bulk fluid, this is not a requirement. Here we show that a generalization of the Ewald approach,^13,14 utilizing a grand canonical treatment with single ion insertions and deletions, allows us to maintain a given surface potential relative to a fixed reference, which will generally not be that of the bulk fluid. Single ion insertion and deletions are not usual in grand canonical schemes for charged fluids, due to obvious problems associated with non-electroneutrality. However, this is avoided in the presence of conducting surfaces by the counterbalancing charge of the images.

A specific example, where we expect our method to be particularly useful, is the theoretical study of electric double layer capacitors, sometimes called “supercapacitors”.^15–19 This area has attracted substantial experimental and theoretical work.^3,15–31 Some experimental data^20,21 support the existence of an electrode pore size (typically of the order of nanometers), which maximizes the capacitance of the system. However, controlling and measuring size distributions of nanoporous electrodes is difficult, and the measurements would strongly benefit by complementing theoretical modelling. Mean-field analyses have indeed suggested the presence of “superionic states” in these systems, whereby a very strong electrostatic screening in narrow pores of conducting materials might facilitate dense packing of equally charged species, which in turn can generate a huge capacitance.³² Still, such simple mean-field predictions need validation by more exact treatments, such as simulations. Simulation studies of porous electrodes are notoriously difficult due to the presence of multiple image charges. Dudka et al. have made attempts to circumvent these problems by simulating lattice models of ionic liquids, where interactions are restricted to nearest-neighbours.³¹ While such models certainly can generate interesting insight, it is nevertheless of considerable interest to be able to utilize a continuum space fluid model, in which Coulomb interactions (including surface polarization) are properly treated.

In the next section, we describe the Ewald method we use to treat the ionic fluid between two conducting surfaces. The original derivation of this model was given by Hautman et al.³³ some three decades earlier, but it seems not to have been widely adopted in recent treatments of these systems. The description contained here has some slight differences, in particular, we define a different unit cell which is useful for charged electrodes. Furthermore, as mentioned earlier, we describe the use of a grand canonical ensemble simulation scheme to create charged electrodes with a potential relative to a well-defined reference.

2 Method

Our system consists of an ionic fluid confined between two infinite and flat perfectly conducting surfaces (electrodes), parallel to the (x, y) plane, and positioned at z = ±h/2. The electrode surface is defined by the dielectric discontinuity, which is assumed to be sharp. In simulations we necessarily restrict the system size to a unit cell, C₀, defined by |x| ≤ L_x/2 and |y| ≤ L_y/2, containing a finite number of ions. The simulation box need not be electro-neutral.

Periodic boundary conditions (PBC) in the (x, y) directions give an infinite number of adjacent replicas of C₀ parallel to the confining electrodes (see Fig. 1(a)), which can be treated using a suitable 2-dimensional (2D) Ewald scheme. The polarization of the two conducting electrodes can be described with image charges, resulting in an infinite number of reflections³⁴ in the two subspaces z > h/2 and z < −h/2. These images give rise to a perpendicular electrostatic field at the dielectric discontinuity, which mimics the equilibrium response of electrons in the conductors. The electrodes thus have a constant potential over their surfaces. The combination of a 2D Ewald scheme with a very large number of explicit image reflections, amounts to a significant computational task for most reasonably sized systems and alternatives to such a straightforward approach are desirable. For example, Santos et al.⁴ have recently described a more efficient but complex treatment for this system, which utilizes periodic Greens functions and a modified 3D Ewald summation, combined with a separate evaluation of the contribution from surface polarizations.

Fig. 1 (a) A unit cell C₀ (b) C₀ reflected across the z = 0 plane. (c) The new unit cell C replicated in the z dimension. Thin arrows indicate periodic boundary conditions while blue and red filled circles represent ions with opposite charge.

We now consider the basic isolated unit cell C₀ (without images). The images are reintroduced by defining a complementary unit cell, C₀′, which consists of the mirror image of C₀ reflected in say the left-hand electrode surface. This process not only reflects the particle positions but also reverses the sign of the charges, so the total charge of both cells always add to zero. A similar complementary cell can be obtained by reflection in the right hand surface of C₀, see Fig. 1(b). We now define a symmetric electroneutral super-cell, C, consisting of the central cell C₀ and the nearest halves of the adjacent complementary cells, see Fig. 1(c). Thus C contains open boundaries, but possesses internal bounding surfaces which correspond to the positions of the original electrodes. Furthermore, though C contains 2N particles, only N are independent, therefore we can write the energy of particles in C as a function of the co-ordinates of particles in the sub-cell C₀.

An infinite 3-dimensional array of replicas of C reproduces the original ionic fluid, between two conducting electrodes, and with PBCs along the (x, y) directions. Such an array can be easily treated with usual 3D Ewald methods to obtain the energy of the ions in C, which is numerically equal to the sum of the energy of a sub-cell C₀ and a complementary cell, C₀′, i.e.,

where U_C0 is the energy of a configuration r^N of ions in C₀ and U_C0′(R[r^N]) is the corresponding energy of the complementary cell with reflected co-ordinates R[r^N] and oppositely signed charges. Due to symmetry, U_C0(r^N) = U_C0′(R[r^N]) and we therefore obtain the simple result for the required energy, U_C0,

We shall denote this approach as the “Image Ewald” method (IE). We note that in the direction parallel to the surfaces the periodicity may induce crystal artifacts, yet by utilizing isotropic periodic boundary conditions³⁵ in these directions such effects will be mitigated. A direct effect from such a scheme is that the reciprocal-space Ewald forces in the xy-directions will be zero, yet there will be a long-range contribution to the total energy. In our simulations, however, we did not use isotropic boundary conditions.

The total Ewald energy for the cell C will have the general form

U_C(r^N) = U_C′(r^N) + U_surf(r^N) (3)

where U_C′(r^N) is a shape-independent part and U_surf(r^N) is the dipolar contribution that gives rise to a shape-dependent surface term. The shape-independent term will give rise to a potential of exactly zero at each electrode, due to the symmetry created by the image charges. When the cells are summed over a spherical geometry, for example, the surface term can be eliminated via tin-foil boundary conditions, which is appropriate in a system of neutral electrodes, that has now been mapped onto a properly three-dimensional lattice.

However, for the general case of charged electrodes, one needs to take cognizance of the slab geometry. This is because, the surface term will be responsible for creating a potential difference between electrodes. For the slab geometry the surface term is given by,

U_surf(r^N) = 2πM_z²/V_cell (4)

where M_z is the total instantaneous dipole moment in the z-direction (perpendicular to the electrodes) and V_cell is the cell volume. This result applies irrespective of the dielectric constant of the medium in which the system is imbedded. That is, the surface term remains even in the presence of tin-foil boundary conditions. We should note, however, that the presence of this uniform polarization field does not violate the condition that each electrode is an equipotential surface. In the case where the charged electrodes are identical, symmetry dictates that C has zero average dipole moment 〈M_z〉 = 0. For a large enough cell, we expect the majority of configurations to have M_z ≈ 0 and therefore it makes physical sense to ignore the surface term in our calculations, as in the case of neutral electrodes. These will be the types of systems considered in this study. As an aside, we also note that the use of 3D Ewald with an inserted slab of vacuum in order to simulate 2D boundary conditions, does require the presence of this surface energy in order to eliminate unwanted electrostatic contributions from neighboring fluid layers. This is not required in the present case where adjacent layers are required to interact. In the case of asymmetric electrodes we will have that, 〈M_z〉 ≠ 0, and we will need to properly account for the potential difference between electrodes that the averaged surface term will bring about.

The IE method described here offers significant computational advantages over many other recent approaches that have been used to deal with such systems, both in terms of efficiency and simplicity of coding. As mentioned earlier, Hautman et al.³³ originally described this approach some years ago, however, that work seems not to have gained wide recognition. Since we developed the methodology without prior knowledge of the earlier work by Hautman et al.,³³ there are actually some conceptual differences. For example, Hautman et al. used a different type of unit cell consisting of C₀ and its adjacent complement C₀′. It is clear that even for symmetric, charged electrodes the average dipole moment in such a cell is not zero, giving a non-zero average surface energy in the slab geometry. This would lead to a spurious electric field within the cell that would need to be removed in an ad hoc fashion. We therefore, argue that our choice of unit cell is conceptually preferable. Furthermore, as far as we are aware, the method of Hautman et al. has not been applied to the case of electrodes with a net total surface charge. This is likely because the original formulation uses a Hamiltonian wherein the electrode potentials are specified. Such an approach is precluded when these potentials (relative to a fixed reference) are unknown. We will now show how the IE can be used in a grand canonical ensemble simulation, which naturally gives rise to net charged electrodes and solves the problem of specifying a fixed reference potential.

2.1 Grand canonical simulations

To simulate ionic fluids bound by electrodes we used the Grand Canonical (Metropolis) Monte Carlo (GCMC) method. When employing the GCMC method it is usual to choose an initially electroneutral configuration with subsequent insertions and deletions of electroneutral combinations of ions. These ions are imagined to be taken to or from a fictitious bulk, and it is useful to define an “average” (experimentally measurable) chemical potential per ion, μ_salt. This can be obtained by averaging the individual chemical potentials, μ_ind, i.e., , where N_sp is the number of species with different valencies Z_i. The coefficients ν_i satisfy and . Individual ion chemical potentials are experimentally elusive to measure, but it is relevant to note that, given the definition of μ_salt, it is possible to write the individual ion chemical potentials as,

μ_ind(Z_i) = μ_salt − Z_ieΨ_ref (5)

where Ψ_ref measures the asymmetry of the model. For a Z : Z electrolyte in the restricted primitive model (RPM) Ψ_ref = 0, (see ESI†) but it is non-zero in a 2 : 1 model, for example. Hence, the IE provides a method to determine μ_ind(Z_i), as we will explicitly demonstrate below, where we also compare IE results with those from the established modified Widom technique.³⁶

In our simulation method, electroneutrality is always maintained by the correlated presence of oppositely charged ions in the complementary C₀′ sub-cell. This is physically equivalent to the creation of counter-balancing image charges in the conducting electrodes and allows us to perform simulations at constant electrode potential, rather than constant electrode charge, using single ion insertions and deletions. Constant potential systems are arguably more relevant to real experimental scenarios. The fluctuating electrode surface charge density, σ_S, will have an average value 〈σ_S〉 that can be regulated by application of a constant potential that will bias either anions or cations via the chemical potentials. With a μ_salt for a given bulk electrolyte, “effective” individual ion chemical potentials are constructed using

μ_eff(Z_i) = μ_salt − Z_ieΨ_bias (6)

where the biasing potential is denoted Ψ_bias. Eqn (6) enters the grand canonical Boltzmann weights for the insertion/deletion moves for the individual ions.

2.2 Simulation details

The ionic fluid was modelled using the RPM wherein the solvent only enters via the dielectric constant, ε_r = 78.3. At a temperature of 298 K, this gives a Bjerrum length of about 7.16 Å. In order to highlight some interesting aspects of asymmetric electrolytes, we assumed a 2 : 1 salt, consisting of divalent and opposing monovalent charges embedded in hard spheres of diameter 5 Å. Generalization to other kinds of salt is straightforward. The Ewald sums used a minimum image cutoff in the real space sum. We set L_x = L_y = 200 Å, and used an equal number of reciprocal wave vectors in the Fourier sum along the reciprocal z as in the reciprocal (x, y) directions, with a suitable symmetrization in the latter case.³⁷ A total of 1800 vectors were used for the Fourier component. The Ewald splitting parameter, often denoted as α, was set to 7/L_x. At least 2 × 10⁷ configurations were used for the final production runs, subsequent to equilibrations, which used even more configurations.

For all simulation results presented in the main article, we have set βμ_salt = −10.7 where β is the inverse thermal energy. In the ESI,† we provide results for other scenarios, including a 3 : 1 salt, and a monovalent salt where the cations carry a non-central charge.

3 Results

A quantity of interest is the average electrostatic potential, Ψ(z), which can be evaluated by integration of the average charge density profile, as described in the ESI.† This calculation simply adds the direct contributions from the charge densities and (due to symmetry) will always give a value of zero, for the potential at the electrode surfaces. Fig. 2(a) and (b) shows the average electrostatic potential Ψ(z) profiles from our GCMC simulations, using IE, with a slit of width, h = 145 Å. This is obviously large enough to produce a bulk fluid in the mid-plane region as evidenced by the flat profile there. Two different Ψ_bias are chosen, and these are compared with the result for Ψ_bias = 0. Normally, we measure Ψ(z) relative to the bulk, which requires addition of a constant potential correction, called the “Donnan potential”,³⁸Ψ_Donn, due to charge distributions external to the slit. Due to the asymmetry of the electrolyte model, Ψ_Donn is not simply equal to Ψ_bias, since there is a finite Ψ_ref in the bulk. From eqn (6) and (5) can write,Therefore we have, Ψ_Donn = Ψ_bias − Ψ_ref. We note that Ψ_Donn is then also equal to the surface potential of the electrodes relative to the bulk fluid. But while Ψ_bias is chosen a priori, Ψ_ref remains undetermined. It can be obtained, however, by noting that the mid-plane potential behaves as,

Ψ(z = 0) → −Ψ_Donn (8)

as h→∞. That is, the potential (relative to the bulk) at the center of the slit should be zero for wide slits as one has effectively a bulk fluid there. From Ψ(z) in Fig. 2 for the case Ψ_bias = 0, we obtain a mid-plane value essentially equal to Ψ_ref, which follows from eqn (8). We then verify that the mid-plane values for the other Ψ(z) do approach −Ψ_Donn.

Fig. 2 Average potential profiles, as obtained from Image Ewald simulations. (a) βΨ_biase = −2. The profile for Ψ_bias = 0 is also shown. (b) βΨ_biase = −3. The profile for Ψ_bias = 0 is also shown.

We evaluated the IE by comparing it with a more straightforward approach, where we use GCMC simulations, as discussed above, but explicitly evaluate a finite number of multiple image reflections between the electrodes. Previous work has noted that convergence of the sum of reflected images is facilitated by using an even number of terms.⁹ Further details are also provided in the ESI.† We also used the minimum image (MI) convention to treat the PBC in the x, y directions, rather than Ewald sums. Hence we use the acronym MI to denote this approach. In the MI method, long-range interactions outside the MI truncation were estimated by an “external” mean-field potential,^39,40 from previously established average charge densities, appropriately symmetrized, and adjusted so that the external potential vanishes at the surfaces. It should be noted that the inclusion of this long-range correction only had a small influence on the properties investigated here, i.e. the results were very similar to those obtained with a simple MI truncation, without any long-range correction. Fig. 3 summarizes our comparisons between the IE and MI approaches for βΨ_biase = −2 and h = 145 Å. As noted above this width is sufficient to produce bulk conditions in the center of the slit.

Fig. 3 Comparisons between Image Ewald results, and MI simulations with a limited number of image reflections, using βΨ_biase = −2. (a) Locally measured mid plane and surface potentials, Ψ_mid, and Ψ_surf. In the MI approach, these are randomly sampled during the simulations. With the Image Ewald method, Ψ_surf = 0, and Ψ_mid is obtained by integration of the average charge density profile. (b) The average surface charge density, 〈σ_S〉. (c) Integral capacitances, C_I.

In Fig. 3(a), we note how the MI method gives a surface value of Ψ(z), which varies significantly from zero, with the number of reflections used. The potential difference between the mid-plane and the surface in the MI does approach the IE result asymptotically but the convergence is quite slow. Surprisingly, these results indicate that the mid-plane potential from the MI actually displays only a modest variation with the number of image reflections, and is in all cases rather close to the value obtained by the IE. One could then be tempted to just assume that Ψ(z) is zero at the surface and use the mid-plane value to set the surface potential relative to the bulk (recall that the negative of the mid-plane value is equal to the proper reference surface potential). However, such a strategy belies other serious errors obtained with too few reflections. For example, in Fig. 3(b), we see that the average surface charge density varies strongly with the number of image reflections, only slowly converging towards the IE value. In Fig. 3(c), we compare predictions of the integral capacitance C_I. Using Image Ewald, this is simply obtained as C_I = 〈Δσ_S〉/Ψ_bias, where 〈Δσ_S〉 is the difference between the average surface charge density at Ψ_bias and the average surface charge density at with a vanishing bias potential (i.e. Ψ_bias = 0). However, with the MI approach, we have to account for truncation errors, which leads to a finite (local) surface potential. Interestingly enough, the capacitance converges faster towards the Image Ewald value than the net potential, or the average surface charge density. This is obviously due to a cancellation of errors. It should be emphasized that since we generally would define the differential capacitance (C_D) as C_D = Δσ_s/ΔΨ_Donn = Δσ_s/ΔΨ_bias, we do not need to establish Ψ_ref, to generate a C_D curve. This is a significant advantage for narrow pores (slits).

We now turn our attention to the individual chemical potentials, μ_ind. Using eqn (5) we calculate μ_ind for the 2 : 1 RPM salt using our method, by setting Ψ_bias = 0 and measuring Ψ_ref. We compare with results from the modified Widom technique using canonical bulk simulations, a cubic box, and Coulomb interactions that are truncated by the minimum image convention. The canonical densities were given by the mid-plane density from grand canonical IE simulations, at zero applied bias potential. Size convergence was estimated by simulating three different simulation box sizes, using the modified Widom method. The results are summarized in Fig. 4, where the standard deviations from the canonical simulations are smaller than the size of the symbols, and therefore not indicated. We note a modest but significant size dependence of the modified Widom simulations, but also that the large size limit seems to approach the Image Ewald result.

Fig. 4 Calculated individual chemical potentials for the cation (top) and anion (bottom). The dashed line indicates the value acquired using Image Ewald simulations, combined with eqn (5).

4 Conclusions

In summary, we have utilized a simple and effective way to manage long-range interactions, using an exact model of infinite image reflections, in systems where an ionic fluid is confined between two charged (in general), planar, and perfectly conducting surfaces, or electrodes. Note that “confined” should be interpreted with some care, since our grand canonical method ensures that the fluid is in chemical equilibrium with a bulk solution, at some prescribed chemical potential. Using wide surface separations, our system will accurately approximate the behaviours at a single electrode surface and we are able to calculate the individual chemical potentials. An important advantage is that we are also able to model a porous electrode, by utilizing a short surface separation. Our method relies on straightforward 3D Ewald techniques, which are wide-spread, and can therefore be highly optimized using particle-mesh Ewald, smooth particle-mesh Ewald etc., to achieve a computational scaling of [scr O, script letter O](N log N), where N is the number of charges in the system. Hence, the Image Ewald method is suitable for a multitude of system setups. One such case is the study of electrochemical behaviours of salt solutions and ionic liquids, in the presence of nanoporous electrodes. Finally, we note that while we have utilized the Metropolis Monte Carlo technique in this work, the IE method can naturally be directly implemented in molecular dynamics codes.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

J. F. acknowledges financial support by the Swedish Research Council.

Notes and references

1. C. Merlet, B. Rotenberg, P. A. Madden, P.-L. Taberna, P. Simon, Y. Gogotsi and M. Salanne, Nat. Mater., 2012, 11, 306–310.
2. E. Frackowiak, Phys. Chem. Chem. Phys., 2007, 9, 1774–1785.
3. S. Kondrat, N. Georgi, M. V. Fedorov and A. A. Kornyshev, Phys. Chem. Chem. Phys., 2011, 13, 11359–11366.
4. A. P. dos Santos, M. Girotto and Y. Levin, J. Chem. Phys., 2017, 147, 184105.
5. S. Tyagi, A. Arnold and C. Holm, J. Chem. Phys., 2007, 127, 154723.
6. S. Tyagi, A. Arnold and C. Holm, J. Chem. Phys., 2008, 129, 204102.
7. R. Kjellander and S. Marcelja, J. Chem. Phys., 1985, 82, 2122–2135.
8. D. Bratko, B. Jönsson and H. Wennerström, Chem. Phys. Lett., 1986, 128, 449–454.
9. R. Szparaga, C. E. Woodward and J. Forsman, Soft Matter, 2015, 11, 4011–4021.
10. S. K. Reed, O. J. Lanning and P. A. Madden, J. Chem. Phys., 2007, 126, 084704.
11. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
12. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
13. P. P. Ewald, Ann. Phys., 1921, 369, 253–287.
14. T. M. Nymand and P. Linse, J. Chem. Phys., 2000, 112, 6152–6160.
15. P. Simon and Y. Gogotsi, Nat. Mater., 2008, 7, 845–854.
16. Y. Bai, Y. Cao, J. Zhang, M. Wang, R. Li, P. Wang, S. M. Zakeeruddin and M. Grätzel, Nat. Mater., 2008, 7, 626–630.
17. S. Ito, S. M. Zakeeruddin, P. Comte, P. Liska, D. Kuang and M. Grätzel, Nat. Photonics, 2008, 2, 693–698.
18. V. Lockett, R. Sedev, J. Ralston, M. Horne and T. Rodopoulos, J. Phys. Chem. C, 2008, 112, 7486–7495.
19. M. Armand, F. Endres, D. R. MacFarlane, H. Ohno and B. Scrosati, Nat. Mater., 2009, 8, 621–629.
20. J. Chimola, G. Yushin, Y. Gogotsi, C. Portret, P. Simon and P. Taberna, Science, 2006, 313, 1760.
21. C. Largeot, C. Portet, J. Chmiola, P.-L. Taberna, Y. Gogotsi and P. Simon, J. Am. Chem. Soc., 2008, 130, 2730–2731.
22. A. A. Kornyshev, J. Phys. Chem. B, 2007, 111, 5545–5557.
23. D. Roy and M. Maroncelli, J. Phys. Chem. B, 2010, 114, 12629–12631.
24. J. Forsman, C. E. Woodward and M. Trulsson, J. Phys. Chem. B, 2011, 115, 4606–4612.
25. D.-e. Jiang, Z. Jin and J. Wu, Nano Lett., 2011, 11, 5373–5377.
26. M. Z. Bazant, B. D. Storey and A. A. Kornyshev, Phys. Rev. Lett., 2011, 106, 046102.
27. C. E. Szparaga, R. Woodward and J. Forsman, J. Phys. Chem. C, 2012, 116, 15946.
28. C. Merlet, M. Salanne and B. Rotenberg, J. Phys. Chem. C, 2012, 116, 7687–7693.
29. C. Merlet, D. T. Limmer, M. Salanne, R. van Roij, P. A. Madden, D. Chandler and B. Rotenberg, J. Phys. Chem. C, 2014, 118, 18291–18298.
30. A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2753–2757.
31. M. Dudka, S. Kondrat, O. Benichou, A. A. Kornyshev and G. Oshanin, J. Chem. Phys., 2019, 151, 184105.
32. S. Kondrat and A. Kornyshev, J. Phys.: Condens. Matter, 2010, 23, 022201.
33. J. Hautman, J. W. Halley and Y. Rhee, J. Chem. Phys., 1989, 91, 467–472.
34. M. M. Hatlo and L. Lue, Soft Matter, 2008, 4, 1582–1596.
35. B. Stenqvist and M. Lund, EPL, 2018, 123, 10003.
36. B. R. Svensson and C. E. Woodward, Mol. Phys., 1988, 64, 247–259.
37. B. A. Wells and A. L. Chaffee, J. Chem. Theory Comput., 2015, 11, 3684–3695.
38. F. G. Donnan, Z. Elektrochem. Angew. Phys. Chem., 1911, 17, 572.
39. G. M. Torrie and J. P. Valleau, J. Chem. Phys., 1980, 73, 5807–5816.
40. B. Jönsson, H. Wennerström and B. Halle, J. Phys. Chem., 1980, 84, 2179.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp01640c

---