A Singlet Reference Interation Site Model Theory ... - ACS Publications

J. Phys. Chem. B 2008, 112, 3365-3374

3365

A Singlet Reference Interation Site Model Theory for Solid/Liquid Interfaces Part II: Electrical Double Layers Stefan Woelki,*,† Hans-Helmut Kohler,† and Hartmut Krienke‡ Institute of Analytical Chemistry, Chemo- and Biosensors, UniVersity of Regensburg, D-93040 Regensburg, Germany, and Institute of Physical and Theoretical Chemistry, UniVersity of Regensburg, D-93040 Regensburg, Germany ReceiVed: September 17, 2007; In Final Form: December 3, 2007

The previously established singlet reference interaction site model (SRISM) theory for the calculation of the fluid structure in the vicinity of a plane impenetrable interface is renormalized for the application to electrical double layers. In combination with the HNC and KH closures, the equations are solved numerically for a 1 M electrolyte solution adjacent to a charged wall with varying surface charge densities. The wall-solvent and wall-ion density distributions as well as the profiles of the electrical field and the electrical potential are compared to computer simulation results. Reasonable agreement is obtained.

1. Introduction The current paper constitutes the second part of a series dealing with a novel integral equation approach to inhomogeneous fluids. In Part I, we established a singlet reference interaction site model (SRISM) theory that allows for the calculation of inhomogeneous density distributions of arbitrary molecular liquids, including mixtures and solutions in the vicinity of a plane, impenetrable wall.1 The theory has been applied to pure water and aqueous electrolyte solutions adjacent to an uncharged soft wall, and the results have shown reasonable agreement with atomistic computer simulations for similar systems. In this paper, the SRISM approach is extended to charged walls, that is, to electrical double layers (EDLs). An overview of numerous theoretical approaches to inhomogeneous fluids including EDLs has already been given in the first part of this series,1 where a number of reviews and continuative references to research papers on this issue can be found. However, irrespective of the individual approach, homogeneous or inhomogeneous fluids involving long-ranged Coulombic forces always make higher demands than systems controlled by purely short-ranged interactions. The pronounced long-range character of Coulombic forces gives rise to a numerically divergent behavior in both integral equation theories and computer simulation methods. To overcome these difficulties, a number of mathematical techniques have been developed based on some appropriate division of the total interaction potential into short-ranged and long-ranged (Coulombic) parts. The latter are treated separately to meet requirements of the individual theoretical or computational method. For computer simulations, a variety of Ewald summation methods have been developed that can be traced back to the pioneering work of Paul Peter Ewald published more than 80 years ago2 and aim to replace the cumbersome summation of electrostatic interaction energies in real space by a rapidly convergent summation in Fourier space.3,4 * To whom correspondence should be addressed. E-mail: stefan.woelki@ chemie.uni-regensburg.de. † Institute of Analytical Chemistry, Chemo- and Biosensors. ‡ Institute of Physical and Theoretical Chemistry.

Similar techniques involving Fourier transforms are also employed to solve Ornstein-Zernicke (OZ)-type integral equation approaches,5 such as RISM6,7 or DRISM8,9 theories, for homogeneous bulk fluids. Such theories are intrinsically tied to and solved in Fourier space, and it is obvious to split and renormalize the total interaction potentials in such a way that the numerically divergent Coulombic parts become amenable to a purely analytical Fourier transformation while the shortranged parts remain well-suited for numerical treatment.10,11 Somewhat different techniques are required to extend singlet integral equation theories for simple fluids in the vicinity of an uncharged wall such as the well-known SOZ theory12-14 (also referred to as OZ1, Henderson-Abraham-Barker or HAB theory) to electrified solid/liquid interfaces. This problem is illsuited for Fourier transform methods.15 Moreover, the unscreened Coulomb potential between a charged wall and an ion is of infinite range, rendering a numerical solution of the SOZ equation in its original notation14 impossible. Henderson and co-workers16 have solved this problem by separating the Coulombic part of the total wall-fluid particle potential under consideration of the global electroneutrality condition. In the following section, we will show that this approach can be adopted for an extension of the SRISM theory to electrified solid/liquid interfaces. The analytical procedures, however, will turn out to be more involved in the case of the SRISM theory than for the SOZ theory. In section 3 we will briefly discuss the solvent and ion models employed in our SRISM calculations. In section 4 the resulting solvent and ion distributions will be presented and compared to corresponding computer simulation results. 2. Renormalization of the SRISM Theory Within the framework of RISM approaches, the structure of a liquid in the vicinity of a solid/liquid interface can be characterized by the normalized density profiles gβ(z) of the fluid’s individual interaction sites β ) 1...K at a distance z from the interface. The functions gβ(z) are defined by the ratio between the local number density Fm β (z) and the average bulk of the fluid molecules containing site(s) β number density Fm β

10.1021/jp077485z CCC: $40.75 © 2008 American Chemical Society Published on Web 02/27/2008

3366 J. Phys. Chem. B, Vol. 112, No. 11, 2008

Woelki et al.

and are related to the singlet direct and indirect wall-fluid site correlation functions, cβ(z) and γβ(z), by

gβ(z) ) γβ(z) + cβ(z) + 1

(2.1)

We have shown in Part I1 that the correlation functions cβ(z) and γβ(z) are related by a symmetry-reduced SRISM equation, which can be written as K

γβ(z) )

∑ νκFmκ ∫-∞ κ)1

+∞

νκ

K

∑ ∑ κ)1 u)1 β

κu*β*

cκ(t)H (βκ)(z - t)dt +

1 2l

β*κu

∫z-lz+l

β*κu

β*κu

cκ(t)dt

β ) 1...K (2.2)

In the notation of eq 2.2, equivalent interaction sites (e.g., the two hydrogen atoms in a water molecule) are combined to a single interaction site group κ with νκ members κ1...κνk and are represented by a single representative site κ*. Such symmetry reduction17,18 avoids the multiple computation of identical correlation functions in the traditional RISM approach. lβ*κu is the intramolecular distance between the representative site β* of group β and the uth site of group κ. The superscript β at the sum in the second term of the right-hand side of eq 2.2 indicates that only groups κ that form part of the same molecule as β are considered. The functions H (βκ) are defined in terms of the bulk site-site total correlation functions h(Rβ):

H (βκ)(V) ) 2π

∫|V|+∞ h(βκ)(s)sds

(2.3)

As usual, eq 2.2 must be solved together with a second type of equation relating γβ(z) and cβ(z). For inhomogeneous solid/liquid systems, such closure relation can generally be written as

cβ(z) )

{

}

exp [-uβ(z)/kT + γβ(z) + Bβ(z)] - γβ(z) - 1 z > 0 z 0 (2.5b) β

β

β

β

In combination with the SRISM equation (eq 2.2), we used these approximations, eqs 2.5a and 2.5b, in Part I of this series to calculate the inhomogeneous wall-fluid site density profiles of pure water and aqueous electrolyte solutions adjacent to an uncharged plane interface.1 Additionally, we employed the modified Duh-Haymet (MDH) approximation,20 which, however, did not yield significant improvements over the HNC and KH closures and, thus, is no longer considered in the present paper. The water sites and the ions were modeled as charged

Lennard-Jones (LJ) spheres with LJ parameters σββ and ββ and were allowed to interact with an infinite monolayer of smearedout uncharged LJ spheres with LJ parameters σ00 and 00, giving rise to the following short-ranged 10 - 4 dispersion potential uβs between the wall and a site β:

uβs (z) ) 4πξ00βσ0β2

[( ) 1 σ0β 5 z

10

-

( )]

1 σ0β 2 z

4

(2.6)

ξ0 is the number density of wall particles per unit area, σ0β ) (σ00 + σββ)/2, and 0β ) (00‚ββ)0.5. From a theoretical point of view, the SRISM approach is in no way restricted to shortranged wall-fluid site potentials (such as the dispersion potential in eq 2.6) but also holds for EDLs involving long-ranged Coulombic interactions uβl between a fluid site β with charge qβ and an electrified interface of charge density σ. In this case, uβl is

uβl ) -

qβσ |z| 2o

(2.7)

o is the vacuum permittivity. The electrostatic potential at the interface (z ) 0) is arbitrarily set to zero. Nevertheless, application of the SRISM theory (eqs 2.2 and 2.4) to EDLs is not straightforward, as is seen from the following argument: Convergence of the individual integrals in the SRISM equation eq 2.2 necessarily requires finite limits limzf-∞ cβ(z) ) cβ- and limzf+∞ cβ(z) ) cβ+. Issues concerning the negative limits cβhave been addressed in detail in Part I of this series1 and will not be discussed here. The positive limits cβ+ can directly be deduced from the closure eq 2.4, which, for large distances z . 0, simplifies to

cβ(z) ≈ -

uβ(z) kT

z.0

(2.8)

For wall-fluid site potentials uβ(z) that asymptotically converge to zero, such as the dispersion potential in eq 2.6, the boundary condition cβ+ ) constant imposes no problems. However, long-ranged wall-fluid site Coulomb potentials ul, given by eq 2.7, do not converge at all, but rather increase or decrease linearly with increasing distance z and, consequently, give rise to strongly divergent terms that render a numerical solution of the SRISM theory in the notation of eq 2.2 with eq 2.4 quite hopeless. To overcome this problem, we will show in the following that the SRISM theory can be rewritten in a renormalized form where the entirety of divergent Coulombic terms completely cancels out as a result of the global electroneutrality condition. First, we separate the divergent long-range Coulombic parts of the direct correlation functions by splitting cβ(z) into a wellbehaved short-range contribution cβs (z) and a long-range part cβl (z):

cβ(z) ) cβs (z) + cβl (z) ) cβs (z) +

qβ σ |z| 2okT

(2.9a)

In a similar manner, we define a short-ranged indirect correlation function γβs (z) by

γβ(z) ) γβs (z) - cβl (z) ) γβs (z) -

qβσ |z| 2okT

(2.9b)

Thus, eq 2.1 and the closure eq 2.4 can be rewritten in terms of

SRISM Theory for Electrical Double Layers

J. Phys. Chem. B, Vol. 112, No. 11, 2008 3367

the well-behaved short-ranged functions cβs (z), uβs (z), and γβs (z):

gβ(z) ) γβs (z) + cβs (z) + 1

(2.10)

cβs (z) )

{

}

exp[-uβs (z)/kT + γβs (z) + Bβ(z)] - γβs (z) - 1 z > 0 -γβs (z) - 1 z 0 p > z and z > 0 p < z and z < 0 p > z and z < 0

}

(2.17)

Thus, using the fact that H (βκ)(p), defined in eq 2.3, is an even -z +∞ H (βκ)(p)dp ) ∫+z H (βκ)(p)dp, the first function, i.e., ∫-∞ integral in eq 2.16, for z > 0 (cases I and II), can be written as

∫ ∫

-

(|z - p| - |z|)H (βκ)(p)dp ) z

-∞ +∞

z

∫ (p)dp + ∫ 2∫

+∞

pH (βκ)(p)dp +

pH (βκ)

(|t| - |z|)dt +

z +∞

z

+∞

∫ (p)dp - 2z ∫

pH (βκ)(p)dp - 2z

pH (βκ)

z

∫-∞

+∞

|z - p|H

(βκ)

(p)dp (2.14)

∫

+∞

-∞

+∞

z +∞

z

(p - z)H (βκ)(p)dp

H (βκ)(p)dp )

H (βκ)(p)dp ) z > 0 (2.18a)

(|z - p| - |z|)H (βκ)(p)dp )

∫ 2z ∫ 2z

z

∫z-l

z+lβ*κu β*κu

K

dt ) νβqβ +

∑ κ)1

β

νκqκ ) Qm β

κ*β

(2.15a) and, thus, turns out to be the net charge Qm β of the molecule

∫ pH (p)dp + ∫ (p)dp + ∫ pH (p)dp + ∫ 2 ∫ (p - |z|)H (p)dp z

H (βκ)(p)dp -

-∞ +∞

|z|

K

κu*β*

1

{

∫

The embraced expression in the first term on the right-hand side of eq 2.14 can be simplified to

qβ +

∑ u)1

qκ

On the other hand, for z < 0 (cases III and IV), we obtain

K

∑

νκ

β

(I) - pH (βκ)(p) (II) (p - 2z)H (βκ)(p) (III) (2z - p)H (βκ)(p) (IV) pH (βκ)(p)

|z|dt ) 0

and substituting the integration variable t by p ) z - t in the second integral of eq 2.13, the functions L β(z) can be rewritten as β

(|z - p| - |z|)H (βκ)(p)dp +

Both integrals in eq 2.16 can further be refined by inspecting the various case differentiations arising from the appearance of moduli in the integrands. With respect to the integrand |z - p| - |z|H (βκ)(z) under the first integral on the right-hand side of eq 2.16, we have the following cases:

|z|dt -

κ)1

β

(2.15b)

(|z - p| - |z|)H (βκ)(p) )

K

K

+∞

K

So far, we have merged all divergent Coulombic terms into a set of apparently ill-behaved functions L β(z) with β ) 1...K. Extending the right-hand side of eq 2.13 by the terms K

νκ F m ∑ κ qκ ∫-∞ κ)1

|t|dt +

(βκ)

H (βκ)(p)dp

Equating eq 2.15a with eq 2.15b and inserting the result into eq 2.14, we obtain

L β(z) )

where the function L β(z) is given by

(s)ds )

K

Furthermore, inserting eqs 2.9a and 2.9b into the SRISM equation eq 2.2 and rearranging the result, we obtain

γβs (z) )

sh

H (βκ)

-∞ ∞

|z| +∞

|z|

+∞

(βκ)

z

(βκ)

+∞

|z|

(βκ)

pH (βκ)(p)dp ) pH (βκ)(p)dp ) z < 0 (2.18b)

Eventually, both cases, z < 0 (eq 2.18b) and z > 0 (eq 2.18a), yield identical well-behaved expressions. Hence, eq 2.18b can be used for both cases. Furthermore, the second integral in eq 2.16 can, after some calculations, explicitely be specified in


Woelki et al.

terms of an analytical function mβ*κu(z), given by

mβ*κu(z) )

∫z-lz+l

β*κu

β*κu

(|t| - |z|)dt )

{

|z| > lβ*κu ) |z| < lβ*κu

0 (|z| - l

β*κu 2

}

(2.19)

so that, with eqs 2.18b and 2.19, the function L β(z) in eq 2.16 can be rewritten as K

L β(z) ) 2

νκFm ∑ κ qκ ∫|z| κ)1

+∞

(p - |z|)H (βκ)(p)dp + νκ

K

∑

κ)1

β

qκ

∑

u)1 κu*β*

mβ*κu(z) 2lβ*κu

(2.20)

In contrast to the original representation of eqs 2.4 and 2.2, the renormalized version of the SRISM theory given by eqs 2.11 and 2.12 with eq 2.20 no longer suffers from divergent Coulombic terms and can readily be applied to EDLs. 3. Models We have solved the renormalized SRISM equation (eq 2.12) in combination with the HNC and KH closures (eq 2.4 with eqs 2.5a and 2.5b, respectively), for a 1 M aqueous electrolyte solution adjacent to a charged plane wall of surface charge density σ ranging from -0.3 C‚m-2 to +0.3 C‚m-2. Each calculation has been carried out on two different levels: First, on the so-called Born-Oppenheimer (BO) level involving explicit RISM models for both ions and water and, second, on a primitive McMillan-Mayer (MM) level, where the water is represented by a dielectric continuum with permittivity ) 78 and the solvent-mediated part of the potential of mean force between the ions is approximated by a dielectrically screened Coulomb potential. Usually, inhomogeneous calculations on the MM level are carried out by means of the well-known SOZ equation.14,16 However, we have used the SRISM equation (eq 2.12), which, in the absence of molecular components, is not identical but equivalent to the SOZ equation. The relation between SRISM and SOZ equations was discussed in detail in Part I of this series.1 For simplicity, we will refer to the results of our MM-level and BO-level SRISM calculations as SOZ and SRISM results, respectively. As input, the SRISM and SOZ calculations require the bulk correlation functions of the 1 M electrolyte solution under investigation. These have been calculated from the symmetryreduced DRISM/KH theory18 on the BO level and from the ordinary OZ/KH5,19 theory on the MM level. In the SRISM calculations, the water molecules are represented by a slightly modified SPC/E water model,22 denoted by SPC/E* in the following. The latter has been proposed by Yoshida et al.23 for use in RISM calculations and complements the interaction between water molecules by some small hydrogen-hydrogen LJ potential completely missing in the original SPC/E water model.22 The cations and anions are modeled by oppositely charged but otherwise identical LJ spheres with LJ parameters being equal to those of the SPC/E (and SPC/E*) oxygen site. The interaction site models for water and the ions and the numerical methods and parameters are described in detail in Part I.1 4. Results and Discussion In this section, we present some results of our SRISM calculations in comparison with the computer simulation experi-

Figure 1. Normalized wall-fluid site density profiles gO(z) (water oxygen) and gH(z) (water hydrogen) in a 1 M aqueous electrolyte solution for increasing positive surface charge densities σ. The circles are the MD simulation results and the solid and dashed lines to the SRISM/KH and SRISM/HNC results, respectively.

ments of Crozier and co-workers,21 who investigated similar but not identical systems by inhomogeneous molecular dynamics (MD) computer simulations. The differences in the systems and their effects on the results were discussed in detail in Part I1 and can be summarized as follows: First, Crozier and coworkers21 employed the original SPC/E water22 as the solvent, while we used the slightly different SPC/E* version.23 Second, the simulation system models the wall by a rectangular grid of charged discrete LJ spheres, while the SRISM approach deals with a smooth, homogeneous interface interacting with the fluid sites via homogeneous short-range 10 - 4 dispersion potentials, given in eq 2.6, and a Coulomb potential of the type presented in eq 2.7 instead of the many-body interaction employed in the simulations. Obviously, these differences will be reflected in the results, and quantitative agreement cannot be expected. Wall-Fluid Site Density Distributions. Figure 1 shows the normalized inhomogeneous density distributions gO(z) and gH(z) of the water oxygen and hydrogen sites obtained from computer simulations (circles) and from the SRISM/KH (solid line) and SRISM/HNC (dashed line) theories for several positive surface charge densities σ. It is seen that the initial peak of the simulation oxygen density profile gO(z) appearing at z ) 2.35 Å for σ ) 0 C‚m-2 is both raised and narrowed with increasing surface charge density σ and slightly shifted toward the interface for approximately 0.1 Å per 0.1 C‚m-2. Such a shift is also

SRISM Theory for Electrical Double Layers observed for the second peak at 5.1 Å (σ ) 0 C m-2) for charge densities σ g 0.2 C‚m-2, but its amplitude remains widely unaffected. The SRISM/KH and SRISM/HNC profiles exhibit the same characteristic features, but the effects are less pronounced. The situation is similar for the hydrogen density profiles gH(z). For σ ) 0 C‚m-2, the hydrogen simulation profiles are characterized by two main peaks at 2.3 and 5.5 Å and two intermediate small shoulders at 3.1 and 4.1 Å. With increasing positive surface charge density, the initial simulation hydrogen peak is raised and narrowed but not shifted. On the other hand, an increase of the surface charge density from 0 to +0.3 C‚m-2 shifts the intermediate shoulders and the second main peak about -0.3 and -0.6 Å, respectively. Furthermore, the simulations exhibit a pronounced amplification of the gH(z) oscillations with increasing surface charge density. Except for the small intermediate peak at 3.1 Å, the simulation results for gH(z) at σ ) 0 C‚m-2 are reasonably reproduced by the SRISM/KH and SRISM/HNC calculations, although the amplitudes of the first peak and intermediate shoulder (3.1 Å) are too small or too high, respectively. The effect of positive surface charge on gH(z) with regard to shift and amplification of the peaks is also observed in the SRISM profiles but, as in the case of the oxygen density profiles, are less pronounced. With regard to the interfacial water structure, these results indicate that, at the uncharged wall, the water molecules in the initial layer are aligned more or less in plane-parallel orientation to the interface. This structure is not shaken up by positive surface charges but is even consolidated with increasing surface charge density. Negative surface charge densities affect the oxygen density profiles in a similar manner as positive surfaces charges. This is shown in the left panel of Figure 2. Again, the gO(z) peaks, as obtained from computer simulations, are both amplified and shifted toward the interface for about 0.3 Å if the surface charge density is decreased from 0 to -0.3 C‚m-2. Again, this behavior is satisfactorily reproduced by both SRISM theories. The most striking effect of surface charge on interfacial water structure is observed in the hydrogen density distribution at negative charge densities (right panel in Figure 2). Both the simulations and the SRISM theories predict the appearance of an additional hydrogen peak already at -0.1 C m-2 located very close to the interface at 1.2 Å, which is completely absent at the uncharged interface. This indicates a complete rearrangement of the water molecules in the initial layer from the plane-parallel alignment at neutral or positively charged interfaces to a perpendicular orientation at negatively charged walls. Both the position and the amplitudes of this charge-induced initial hydrogen peak are well reproduced by the SRISM approach, but the positions of the second SRISM/KH and SRISM/HNC hydrogen peaks at negatively charged interfaces differ by about +0.6 Å from the simulation results. However, in view of the approximate character of the SRISM approach and differences in the water and wall models employed in the simulations and in our calculations, the effects of both positive and negative surface charges on the interfacial water structure, as predicted from the SRISM/KH and SRISM/HNC theories, are in reasonable agreement with the reference simulation results. Interestingly, in the case of water, the KH and HNC closures show only insignificant differences in the resulting density profiles. The situation is somewhat different for the ionic density profiles at charged interfaces. Figure 3 shows the SRISM/KH and SRISM/HNC cationic and anionic density profiles g+(z) and g-(z) at positively charged walls together with the corre-


Figure 2. Same as Figure 1 but for increasing negative surface charge densities σ.

sponding computer simulation results by Crozier et al.21 To begin with the uncharged interface as a reference point (upper left and right panel), it is seen that the theoretical SRISM/KH and SRISM/HNC ion distributions exhibit distinct ionic adsorption peaks at the wall-ion contact distance for both cations (2.30 Å) and anions (2.15 Å). Even under consideration of certain ergodicity problems observed in the simulation profiles, such ionic adsorption peaks are completely missing in the simulation results. At positively charged interfaces, one would intuitively expect some accumulation of anions (counterions), while the cations (coions) should be depleted from the interface. Although such charge separation is clearly observed in both the simulation profiles and the SRISM results, there are interesting differences not only between the simulations and the integral equation approaches but also between the SRISM/KH and SRISM/HNC theories: At a moderate surface charge density of σ ) +0.1 C‚m-2, the simulations show practically no counterion (anion) accumulation at contact distance. In contrast to that, the SRISM/ KH theory predicts a pronounced increase of the direct adsorption peak that is even surpassed by the SRISM/HNC approximation. Apart from differences in the initial peaks, the SRISM and simulation counterion profiles g-(z) show reasonable agreement. With regard to the coion (cations), the simulations exhibit the expected coion depletion, which is well reproduced by the integral equation approach, and even more as the SRISM/ HNC and SRISM/KH direct adsorption peaks, occurring at σ ) 0 C‚m-2, are strongly damped at σ ) +0.1 C‚m-2.


Figure 3. Normalized wall-ion density profiles g+(z) (cation) and g-(z) (anion) in a 1 M aqueous electrolyte solution for increasing positive surface charge densities σ. The circles are the MD simulation results, and the solid and dashed lines are the SRISM/KH and SRISM/ HNC results, respectively.

If the surface charge is increased to +0.2 C‚m-2 and further to +0.3 C m-2, the simulation counterion density profiles g-(z) exhibit more pronounced short-range structures characterized by two distinct peaks at 2.1 and 3.8 Å, indicating the formation of two consecutive counterion (anion) layers covering the region 1.0 Å e z e 5.5 Å adjacent to the interface. The distance between the two counterion peaks of about 1 Å is quite small compared to the ion-ion LJ contact distance of 3.166 Å and indicates a considerable overlap of the ion layers. Both peaks are reproduced by the SRISM theories. although there are small deviations of about (0.2 Å in the peak positions and considerable differences concerning the peak heights, especially in the case of the HNC approximation. In other words, both SRISM theories and the simulation method agree with respect to the formation of a 2-fold counterion layer at the positively charged wall but differ in the resulting partitioning of counterions and coions among the two layers. The SRISM/HNC theory puts a vast amount of counterions (anions) into the first layer that almost completely compensate (+0.2 C‚m-2) or even overcompensate (+0.3 C‚m-2) the surface charge density. As a consequence, the second ionic layer contains not only counterions but also a few (+0.2 C‚m-2) or even a majority of coions (+0.3 C‚m-2) to satisfy the global electroneutrality condition. The extreme SRISM/HNC counterion accumulation together with an associated coion adsorption at σ g 0.2 C‚m-2 is due to

Woelki et al.


a well-known shortcoming of the HNC approximation that regularly overestimates the effect of Coulomb interactions.19 In the simulations, the anions are distributed more uniformly among the two counterion layers. Both layers contain anions exclusively and screen approximately two-thirds of the surface charge densities +0.2 C‚m-2 and +0.3 C‚m-2. Outside the initial counterion layers, the simulations show the typical asymptotic behavior of both counterion and coion density profiles. The SRISM/KH theory yields quite similar results: the amount of counterions in the initial layers and the resulting screening effect is close to the simulations, although the layer structure is somewhat softened. The second counterion layer is considerably broader and opens out into the typical asymptotic behavior. Figure 4 shows the situation for a negatively charged interface where the cations and anions act as counterions and coions, respectively. The comparison with Figure 3 shows that the respective counter- and coion profiles clearly depend on the sign of the surface charge density. This effect is due to the asymmetric behavior of the water molecules at positively and negatively charged interfaces and, as will be seen in the following, does not occur in continuum models of EDLs. However, the differences mainly pertain to the intensity of the characteristic effects, while the basic phenomena, such as the formation of 2-fold counterion layers and the HNC-induced recharging of the interface at high charge densities in combination with coion adsorption, are quite similar. Generally speaking, the agreement between the SRISM ionic density profiles and the corresponding simulation results is

SRISM Theory for Electrical Double Layers acceptable albeit less satisfying than the case of the water oxygen and hydrogen density profiles. There are two major deviations: First, the SRISM/HNC theory exhibits some spurious coion adsorption effect at high surface charge densities |σ| g 0.2 C‚m-2, which is completely missing in the simulation results. This effect is due to a characteristic failure of the HNC approximation for systems involving Coulomb interactions and does not appear in the SRISM/KH results. Second, similar to a recent work by Tanimura and co-workers,24 both SRISM theories exhibit direct counter- and coion adsorption at low surface charge densities (|σ| e 0.1 C‚m-2), while the reference simulations predict a complete depletion of the ions from the interface. As discussed in detail in Part I,1 these different ionic adsorption characteristics at low surface charge densities may have several causes: On the one hand, the SRISM direct adsorption peaks may simply be artefacts resulting from the approximate character of the HNC and KH closure relations. On the other hand, ionic adsorption may indeed be a physically sound phenomenon in the case of the smooth wall model employed in the SRISM approach, while the granular wall model used in the reference simulations may contrariwise cause depletion of the ions (see the introductory remarks to this section). As things are now, this point cannot be finally clarified, and more careful investigations are required to draw definite conclusions. Molecular Solvent vs Dielectric Continuum. To estimate the effect of the molecular water model, we have recalculated the ionic density profiles on the MM level (SOZ calculations), where the explicit water molecules are replaced by a dielectric continuum with ) 78. Figure 5 shows the SRISM/KH (left panel) and SRISM/HNC (right panel) counterion density profiles for cations at negatively charged walls (dashed lines), anions at positively charged walls (dashed-dotted lines) and, for comparison, the corresponding SOZ counterion distributions (solid lines). As both cations and anion differ only in their charge, the SOZ counterion profiles are totally symmetric with respect to the sign of the surface charge density and represent both the cation/anode and anion/ cathode cases. It is seen that, in the case of the KH approximation (left panel), the SOZ/KH counterion densities are of the same order of magnitude as the SRISM/KH ion distributions throughout. Generally, the SRISM/KH counterion profiles oscillate around the SOZ/KH curve that, in contrast, is totally smooth and does not exhibit any oscillations. Thus, the formation of counterion (and coion) layers in the immediate vicinity of an interfaces, seen in the SRISM and simulation results (Figures 3 and 4), are obviously due to ion-solvent correlations. At an uncharged wall, the direct ionic adsorption turns out to be more pronounced in the SOZ continuum solvent model than in the explicit RISM water model employed in the SRISM approach. With increasing surface charge density, however, the initial SRISM/KH counterion peaks become increasingly higher than the maximum of the SOZ/KH profiles. The situation is quite similar in the case of the HNC approximation (right panel) but the differences are more pronounced, not only between the SRISM/HNC and SOZ/HNC theories but also with respect to the SRISM anion/cathode and cation/anode asymmetry. Note that the apparent bend in the SOZ/HNC profiles at σ ) (0.2 and (0.3 C‚m-2 is not real but is due to the horizontally split diagrams. With respect to the solvent-induced oscillatory behavior of the SRISM counterion density profiles, it is worth mentioning that such layering effects have not only been observed in


Figure 5. SRISM/KH (left panel) and SRISM/HNC (right panel) counterion density profiles g((z) in a 1 M aqueous electrolyte solution for increasing surface charge densities σ. The dashed and dasheddotted lines are the cation and anion density profiles for negatively and positively charged walls, respectively. The solid lines are the corresponding SOZ/KH (left panel) and SOZ/HNC (right panel) counterion density profiles.

alternative double layer theories involving dipolar25-35 or multipolar36-39 hard sphere water models or the central force water model40-43 as solvent, but surprisingly occur in macroscopic-thermodynamic continuum models as well if both the finite volumes of ions and water and the dielectric saturation of the latter are taken into account.44-46 The corresponding SRISM/KH (left panel) and SRISM/HNC (right panel) coion density profiles for cations at cathodes (dashed lines), anions at anodes (dashed-dotted lines), and the corresponding SOZ results (solid lines) are shown in Figure 6. As in the case of the counterion profiles in Figure 5, the SOZ/ KH and SOZ/HNC coion profiles are smooth and, especially at σ g (0.2 C‚m-2, seem to differ from the SRISM/HNC and SRISM/KH profiles only by the lack of superposed oscillations. Interestingly, at charged interfaces with σ g 0.2 C‚m-2, the SRISM cations/cathode and anions/anode oscillations show a maximum phase shift. Even in the extreme case of σ ) -0.3 C‚m-2, the anions are not completely displaced from the wall and always exhibit a small but visible peak at the wall-anion contact distance (2.3Å). In contrast, the initial cation/cathode peak is considerably shifted to approximately 3.5 Å if the surface charge is increased from 0 to +0.1 C‚m-2. We cannot offer a sound physical explanation for this effect, but, apparently, the integration of anions into an interfacial water layer is entropically


Figure 6. SRISM/KH (left panel) and SRISM/HNC (right panel) coion density profiles g((z) in a 1 M aqueous electrolyte solution for increasing surface charge densities σ. The dashed and dashed-dotted lines are the cation and anion density profiles for positively and negatively charged walls, respectively. The solid lines are the corresponding SOZ/KH (left panel) and SOZ/HNC results.

favored compared to that of cations. However, with regard to that, the SOZ profiles are more similar to the SRISM anion/ anode picture. For both counter- and coions, the structural differences in the ionic density profiles between the SRISM theories and the SOZ continuum approach widely disappear above approximately10 Å. Integral Double Layer Characteristics. Figures 7 and 8 show the SRISM/KH (solid lines) and SRISM/HNC (dashed lines) electric field and potential profiles, E(z) and φ(z) as calculated from the singlet water and ion correlation functions gO(z), gH(z), g+(z), and g-(z). The corresponding formulas are given in Part I.1 For comparison, we have also included the corresponding simulation results (circles) from Crozier et al.21 As to positive surface charges, it is seen in Figure 7 that the integral equation approaches yield quite reasonable agreement with the simulation results. To our surprise, the accordance gets even better with increasing surface charge density. Despite the considerable differences between the SRISM/HNC and SRISM/ KH counterion and coion density profiles, the resulting SRISM/ KH and SRISM/HNC field and potential profiles are quite similar. This is due to fact that, first, the extreme HNC counterion adsorption is always balanced by a likewise accumulation of coions so that the HNC and KH net ionic space charge densities in the fluid are not as different as implied by

Woelki et al.

Figure 7. Electrostatic field profiles E(z) (left panel) and electrostatic potential profiles φ(z) (right panel) in a 1 M aqueous electrolyte solution for increasing positive surface charge densities σ. The circles are the MD simulation results, and the solid and dashed lines are the SRISM/ KH and SRISM/HNC results, respectively.

the individual ion density profiles. Second, because of the overwhelming surplus of water compared to the ions, the total space charge density profiles are anyway much more determined by the oxygen and hydrogen distributions than by the ionic species. With regard to negatively charged interfaces (Figure 8), the deviations of the SRISM/KH and SRISM/HNC electric field and potential profiles from the simulation results are slightly more pronounced than in the case of positively charged walls but are still quite acceptable. 5. Conclusions The current paper deals with the renormalization of the general symmetry-reduced SRISM equation (eq 2.2) derived in Part I of this series1 with respect to Coulombic wall-fluid site interactions. It is shown that the SRISM renormalization can be carried out on an exact analytical level in a way similar to the renormalization of the SOZ equation16 for the treatment of EDLs on the MM level. We have solved the resulting EDL-adapted SRISM equation (eq 2.12) numerically in combination with both the HNC and the KH closure (eqs 2.5a and 2.5b), respectively, for the same 1M model electrolyte solution investigated in Part I1 but, this time, for variously charged interfaces. Even at high surface charge densities, convergence was rapidly obtained for both the

SRISM Theory for Electrical Double Layers

J. Phys. Chem. B, Vol. 112, No. 11, 2008 3373 model40-43 as solvent molecules, the SRISM method has the distinct advantage of much higher flexibility with regard to the chemical components forming the fluid phase. It can readily be applied to EDLs containing mixed solvents, such as water/ acetonitrile or water/acetone mixtures, and molecular ions, such as nitrate and sulfate. Moreover, the SRISM approach is not restricted to the case of one-side confined fluids treated in the current paper, but, following the basic geometric arguments of Lozada-Cassou,47 can also be extended to fluids at nanocylinders or confined in slit pores and nano tubes. Acknowledgment. The authors wish to express their sincere thanks to Prof. Douglas Henderson, Prof. Andriy Kovalenko, and Prof. Georg Schmeer for valuable discussions. The simulation data provided by Dr. Paul Crozier, Dr. Richard Rowley, and Prof. Henderson and control calculations provided by Prof. Georg Schmeer are also gratefully acknowledged. References and Notes


HNC and the KH approximation from initial estimates γβ(z) ) 0. In every case, the computational times were less than 5 min on a standard personal computer, which, in part, is due to the symmetry-reduced notation of the SRISM equation (eq 2.12). The resulting wall-water and wall-ion density distributions as well as the field and potential profiles have been compared to corresponding MD computer simulation results for similar, but not identical systems. Both the SRISM/HNC and the SRISM/KH theories predict surface charge-dependent interfacial phenomena, such as a reordering of the water molecules in the initial water layer at moderate negative surface charges and the formation of 2-fold counterion layers at charged interfaces, which are also observed in the simulation results. Under consideration of the approximate character of the integral equation approach and the differences in the wall and water models employed in the simulations and in the SRISM theories (see the introductory remarks in section 4), the overall agreement is quite acceptable. However, with respect to the resulting counterion distributions at charged interfaces, the KH closure yields better results than the HNC approximation, as the latter widely overestimates the counterion accumulation in the immediate vicinity of the wall. In summary, we think that the EDL-adapted SRISM integral equation theory provides an effective and reliable tool for the investigation of EDLs on a detailed molecular level. Compared to other theoretical BO-level approaches to EDLs, such as inhomogeneous MOZ or OZ theories employing dipolar25-35 or multipolar36-39 hard spheres or the central-force water

(1) Woelki, S.; Kohler, H.-H.; Krienke, H. J. Phys. Chem. B 2007, 111 (47), 13386. (2) Ewald, P. Ann. Phys. 1921, 64, 253. (3) Allen M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1987. (4) Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide; Springer-Verlag: New York, 2002 ;Vol 21. (5) Ornstein, L. S.; Zernike, F. Proc. Acad. Sci. 1914, 17, 793. (6) Chandler, D.; Anderson, H. C. J. Chem. Phys. 1972, 57, 1930. (7) Chandler, D. J. Chem. Phys. 1973, 59, 2742. (8) Perkyns, J. S.; Pettit, B. M. Chem. Phys. Lett. 1992, 190, 626. (9) Perkyns, J. S.; Pettit, B. M. J. Chem. Phys. 1992, 97, 7656. (10) Ng, K. J. Chem. Phys. 1974, 61, 2680. (11) Kovalenko, S.; Ten-No, Hirata, F. J. Comput. Chem. 1999, 20, 928. (12) Percus, J. K. J. Stat. Phys. 1976, 15, 423. (13) Perram, J. W.; Smith, E. R. Chem. Phys. Lett. 1976, 39, 328. (14) Henderson, D.; Abraham, F. F.; Barker, J. A. Mol. Phys. 1976, 31, 1291. (15) Attard, P.; Berard, D. R.; Ursenbach, C. P.; Patey, G. N. Phys. ReV. A 1991, 44, 8224. (16) Henderson, D.; Blum, L.; Smith, W. R. Chem. Phys. Lett. 1979, 63, 381. (17) Bertagnolli, H.; Hausleithner, I.; Steinhauser, O. Chem. Phys. Lett. 1985, 116, 465. (18) Woelki, S.; Kohler, H.-H.; Krienke, H.; Schmeer, G. Phys. Chem. Chem. Phys. 2008, 10, 898. (19) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 1009. (20) Duh, D. M.; Henderson, D. J. Chem. Phys. 1996, 104, 6742. (21) Crozier, P. S.; Rowley, R. L.; Henderson, D. J. Chem. Phys. 2000, 113, 9202. (22) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (23) Yoshida, K.; Yamaguchi, T.; Kovalenko, A.; Hirata, F. J. Phys. Chem. B 2002, 106, 5042. (24) Tanimura, A.; Kovalenko, A.; Hirata, F. Langmuir 2007, 23, 1507. (25) Carnie, S. L.; Chan, D. Y. C. J. Chem. Phys. 1980, 73, 2949. (26) Eggebrecht, J. M.; Isbister, D. J.; Rasaiah, J. C. J. Chem. Phys. 1980, 73, 3980. (27) Henderson, D.; Blum, L. J. Chem. Phys. 1981, 74, 1902. (28) Rasaiah, J. C.; Ibister, D. J.; Stell, G. J. Chem. Phys. 1981, 75, 4707. (29) Dong, W.; Rosinberg, M. L.; Perera, A.; Patey, G. N. J. Chem. Phys. 1988, 89, 4994. (30) Torrie, G. M.; Perera, A.; Patey, G. N. Mol. Phys. 1989, 67, 1337. (31) Berard, D. R.; Patey, G. N. J. Chem. Phys. 1992, 97, 4372. (32) Berard, D. R.; Kinoshita, M.; Ye, X.; Patey, G. N. J. Chem. Phys. 1994, 101, 6271. (33) Berard, D. R.; Kinoshita, M.; Ye, X.; Patey, G. N. J. Chem. Phys. 1995, 102, 1024. (34) Berard, D. R.; Kinoshita, M.; Cann, N. M.; Patey, G. N. J. Chem. Phys. 1997, 107, 4719. (35) Yamamoto, M.; Kinoshita, M. Chem. Phys. Lett. 1997, 274, 513. (36) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1988, 88, 7826. (37) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1988, 89, 3285. (38) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1989, 90, 4513.

3374 J. Phys. Chem. B, Vol. 112, No. 11, 2008 (39) (40) (41) (42) 2792. (43) 94, 1.

Patey, G. N.; Torrie, G. M. Chem. Scr. 1989, 29A, 39. Vossen, M.; Forstmann, F. J. Chem. Phys. 1994, 101, 2379. Vossen, M.; Forstmann, F. Mol. Phys. 1995, 86, 1493. Krmer, A.; Vossen, M.; Forstmann, F. J. Chem. Phys. 1997, 106, Vossen, M.; Forstmann, F.; Kra¨mer, A. Solid State Ionics 1997,

Woelki et al. (44) Woelki, S.; Kohler, H.-H. Chem. Phys. 2000, 261, 411. (45) Woelki, S.; Kohler, H.-H. Chem. Phys. 2000, 261, 421. (46) Woelki, Kohler, H.-H. Chem. Phys. 2004, 306, 209. (47) Louzada-Cassou, M. In Fundamentals of Inhomogenous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; pp 303-361.

A Singlet Reference Interation Site Model Theory ... - ACS Publications

Recommend Documents