Samuel
Stenberg
a,
Björn
Stenqvist
b,
Cliff
Woodward
c and
Jan
Forsman
*a
aTheoretical Chemistry, Naturvetarvägen 14, 22362 Lund, Sweden. E-mail: samuel.stenberg@teokem.lu.se; jan.forsman@teokem.lu.se
bPhysical Chemistry, Naturvetarvägen 14, 22362 Lund, Sweden. E-mail: bjorn.stenqvist@teokem.lu.se
cUniversity College, University of New South Wales (ADFA), Canberra ACT 2600, Australia. E-mail: c.woodward@adfa.edu.au
First published on 2nd June 2020
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.
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 data20,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.32 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.31 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.33 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.
Periodic boundary conditions (PBC) in the (x, y) directions give an infinite number of adjacent replicas of C0 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 reflections34 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.4 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.
We now consider the basic isolated unit cell C0 (without images). The images are reintroduced by defining a complementary unit cell, C0′, which consists of the mirror image of C0 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 C0, see Fig. 1(b). We now define a symmetric electroneutral super-cell, C, consisting of the central cell C0 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 C0.
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 C0 and a complementary cell, C0′, i.e.,
(1) |
(2) |
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 conditions35 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
UC(rN) = UC′(rN) + Usurf(rN) | (3) |
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,
Usurf(rN) = 2πMz2/Vcell | (4) |
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.33 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.,33 there are actually some conceptual differences. For example, Hautman et al. used a different type of unit cell consisting of C0 and its adjacent complement C0′. 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.
μind(Zi) = μsalt − ZieΨref | (5) |
In our simulation method, electroneutrality is always maintained by the correlated presence of oppositely charged ions in the complementary C0′ 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(Zi) = μsalt − ZieΨbias | (6) |
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.
(7) |
Ψ(z = 0) → −ΨDonn | (8) |
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.9 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.
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 CI. Using Image Ewald, this is simply obtained as CI = 〈Δσ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 (CD) as CD = Δσs/ΔΨDonn = Δσs/ΔΨbias, we do not need to establish Ψref, to generate a CD 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). |
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp01640c |
This journal is © the Owner Societies 2020 |