Liquid Interfaces Part I: Uncharged

13386

J. Phys. Chem. B 2007, 111, 13386-13397

A Singlet-RISM Theory for Solid/Liquid Interfaces Part I: Uncharged Walls 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: December 28, 2006; In Final Form: August 23, 2007

We present a novel integral equation method for the calculation of fluid structure in the vicinity of a plane impenetrable wall. The theory is based on the well-known RISM equation and is capable of dealing with arbitrary interaction site model (ISM) fluids at a solid/liquid interface. In conjunction with several closure approximations, the equations are solved numerically and wall-fluid site density distributions as well as charge density, field, and potential profiles are calculated for pure water and aqueous electrolyte solutions with varying concentrations adjacent to an uncharged soft wall. The results show reasonable agreement with corresponding computer simulation data.

1. Introduction Almost 100 years ago, Gouy1 and Chapman2 established the first quantitative theory for electrical double layers (EDLs) on the basis of a macroscopic-thermodynamic Poisson-Boltzmann (PB) approach.3 Since then, a lot of work has been done not only on improvements of the classical Gouy-Chapman model by a number of extended and modified Poisson-Boltzmann (MPB) theories (see refs 4-6 where a large number of further references on MPB theories can be found) but also on the development of more fundamental statistical-mechanical techniques dealing with EDLs and solid/liquid interfaces in general. For details on the various theoretical approaches, the reader is referred to a number of excellent reviews.7-11 The focus of the overwhelming majority of these theories is on so-called primitive models of the liquid phase, where the solvent is described as a dielectric continuum and effects arising from the discrete nature of the solvent particles are completely ignored. However, a number of computer simulations12-21 and experimental investigations, such as surface force measurements,22-26 have shown that the structure of a solvent at an interface differs strongly from that of the bulk phase and, consequently, affects the spatial density profiles of the solutes. In early non-primitive theories for solid/liquid interfaces, the molecular Ornstein-Zernike (MOZ) equation was solved for mixtures of dipolar and charged hard spheres at a plane, charged wall.27-37 The MOZ approach has been further extended to more sophisticated water models by the work of Patey et al.38-41 who investigated mixtures of multipolar and charged hard spheres in the vicinity of charged macrospheres. More recently, Vossen et al.42-45 have employed the central-force water model46 to solve the inhomogeneous Ornstein-Zernike equation (IOZ) for pure water and aqueous electrolyte solutions at a plane charged interface. These calculations have provided important insight into the crucial role of the solvent in EDLs. However, they are restricted to aqueous solutions. For EDLs containing more complex molecular components, the MOZ as well as the CF * To whom correspondence should be addressed. E-mail: stefan. [email protected]. † Institute of Analytical Chemistry, Chemo- and Biosensors. ‡ Institute of Physical and Theoretical Chemistry.

approach becomes extremely cumbersome, not only for the solid/liquid interface but also for the bulk phase. As an alternative to MOZ or CF theories, the reference interaction site model (RISM) theory provides a more widely applicable alternative approach for the description of molecular fluids. In their pioneering work, Chandler and Andersen47,48 have derived a RISM integral equation that, in conjunction with some closure relation, allows for the calculation of orientationaveraged site-site correlation functions between individual interation sites of molecular species. Later, a generalized threedimensional version (3DRISM) has been established.49,50 In its basic version, the RISM theory fails for ionic solutions with finite salt concentrations. It has been shown that this failure is due to inherent dielectric inconsistencies,51,52 and Perkyns and Pettit introduced a correction leading to the dielectrically consistent RISM theory51,53 (DRISM-theory). Hirata has given an extensive review of the RISM approach.54 Although the RISM approach in its various flavors has proven quite successful in a number of investigations on chemically and biologically relevant bulk systems,55-63 there have been remarkably few attempts to extend the RISM approach to inhomogeneous systems and EDLs. Earlier attempts investigate the structure of dumbbells,64 sulfur dioxide,65 and water66,67 near a plane wall in terms of a singlet RISM theory, where the wall is considered as an infinitely diluted, mono-atomic component whose radius approaches infinity. This basic concept of wall implementation has been developed independently by Percus68 and Perram et al.69 for uncharged walls, and in a more general manner, for EDLs on a primitive level in the pioneering work of Henderson et al.70 who established the well-known singlet Ornstein-Zernike (SOZ) equation for primitive isolated EDLs. However, in contrast to the SOZ theory, the limiting process of letting the radius of the wall particle approach infinity has never been carried out analytically in existing singlet RISM approaches.64-67 Instead, the wall-site density distributions were calculated for a series of increasing wall particle radii and the results extrapolated to plane walls. Because of numerical instabilities at large radii, this approach was not very successful. More recently, Akiyama and Hirata71 employed the RISM approach to calculate radial distribution functions between fixed-

10.1021/jp068998t CCC: $37.00 © 2007 American Chemical Society Published on Web 11/06/2007

Singlet-RISM Theory for Solid/Liquid Interfaces in-space surface atoms of a finite sized uncharged crystalline slab and the interaction sites of water. To obtain the corresponding density profiles perpendicular to the interface, Kovalenko and co-workers have solved the 3DRISM theory for water at similiar cristalline72,73 and metal74,75 slabs. The present paper deals with a simpler and more intuitive route for the calculation of density distributions of arbitrary interaction site model fluids, including mixtures and solutions, in the vicinity of a plane wall: We transform the ordinary RISM equation, describing bulk site-site correlations of molecular fluids, into a novel singlet RISM (SRISM) equation that directly addresses the total and direct singlet correlations between a plane wall and the interaction sites of a fluid. In contrast to existing RISM approaches for inhomogenous fluids,64-67 the RISM f SRISM transformation, including the limiting process of inflating a spherical particle to infinite size, is carried out on an exact analytical level avoiding expensive repeated numerical solutions of the RISM equation and uncertain extrapolations to the limit of a plane wall. Compared to inhomogenous MOZ or OZ theories where the solvent molecules are described as dipolar hard spheres or in terms of central-force models, respectively, the SRISM method is not restricted to relatively simple fluids, such as water, but can readily be applied to more complex inhomogenous solid/liquid systems including organic (co-)solvents and solutes. This paper is arranged as follows: In section 2, we show how the RISM equation for site-site correlations in the bulk phase can analytically be transformed to a singlet RISM (SRISM) equation for wall-fluid site correlations. In section 3, we present some useful wall-fluid site interaction potentials and briefly discuss the models employed for the system components in section 4. Section 5 describes the procedure for the numerical solution of the SRISM equations. In section 6, some numerical results for water and electrolyte solutions at interfaces are presented and compared to results from computer simulations for similiar systems.

J. Phys. Chem. B, Vol. 111, No. 47, 2007 13387 group. Note that mono-atomic components are represented as groups with νR ) 1. The orientation averaged radial intermolecular correlation functions between representative reference sites of any two groups R and β are described by the RISM formalism. In this symmetry-reduced component notation, the RISM equation can be written as76 K

hˆ (Rβ) )

∑

K

∑cˆ

ω ˆ (Rλ)

λ)0

(λκ)

(ω ˆ (βκ) + νκFm ˆ (βκ)) κh

κ)0

R ) 0..K, β ) R..K (2.1)

Fm κ is the (bulk) number density of the molecular species containing (νκ equivalent sites of) interaction site group κ. Note that eq 2.1 differs from the usual RISM formalism by the appearance of the factors νκ and a somewhat different definition of the intramolecular correlation functions ω ˆ (Rβ) (see below). In the form of eq 2.1, the redundant calculation of equivalent correlation functions is avoided. For details, the reader is referred to the work of Bertagnolli et al.76 The functions cˆ (Rβ)(k) and hˆ (Rβ)(k) are the Fourier-Bessel transforms of the direct and total radial correlation functions, c(Rβ)(r) and h(Rβ)(r), respectively:

∫c

(r)r sin(kr) dr

∫h

(r)r sin(kr) dr

cˆ (Rβ)(k) )

4π k

hˆ (Rβ)(k) )

4π k

∞ (Rβ)

0

∞ (Rβ)

0

(2.2a) (2.2b)

The total correlation functions h(Rβ)(r) are related to the local density F(Rβ)(r) of the reference site of group β at a distance r from the reference site of R by (Rβ) F(Rβ)(r) ) Fm (r) + 1) β (h

(2.3)

The intramolecular correlation functions ω ˆ (Rβ)(k) appearing in the RISM formalism of eq 2.1 are defined as ω ˆ (Rβ) ) 0

2. Singlet RISM Theory Fundamental Relations. The basic concept of a singlet approach for the description of fluid struture in the vicinity of a plane impenetrable wall is to model the latter as a giant sphere immersed into the bulk fluid. Formally, the wall particle is considered as a component of vanishing number density and infinite radius. To establish a suitable singlet RISM integral equation theory for wall-particle correlations, the RISM equation is modified according to the following recipe: 1. Add an additional spherical component, the “wall” species, to the bulk fluid. 2. Let the number density of the wall species approach zero. 3. Let the radius of the wall species approach infinity. We consider a bulk fluid composed of one or more molecular (or atomic) species. Within the RISM approach, the molecules are constructed as rigid arrangements of interaction sites. According to Bertagnolli et al.,76 the complexity of the system can be reduced by combining all sites of a given molecule that are equivalent under the aspect of symmetry to a single interaction site group. Such symmetry based reduction is of considerable interest as it can lead to important numerical simplifications. Thus, a fluid is formally described by K + 1 distinct interacion site groups R ) 0..K. Group 0 contains only one atom representing the wall component (step 1) while each of the remaining groups R ) 1..K may consist of νR g 1 equivalent interaction sites. In the case of νR > 1, one of these sites is arbitrarily chosen as representative of the respective

R and β do not form part of a common molecule (2.4a) ω ˆ (Rβ) )

νβ

sin klR*βu

u)1

klR*βu

∑

R and β do form part of a common molecule (2.4b) νR

ω ˆ

(RR)

)1+

∑

u)1 Ru*R/

sin klR*Ru klR*Ru R ) β (2.4c)

where lR*βu is the intramolecular distance between the reference site R* of group R and the uth site of group β with 1 e u e νβ. The special case βu ) R* in eq 2.4c is obtained from eq 2.4b ˆ (Rβ) are not by forming the limit lR*βu f 0. Note that the ω ˆ (βR). By defintion, the summation is carried symmetric: ω ˆ (Rβ) * ω out over members of the second index, while the reference site (first index) is fixed. The wall component 0 consists of simple single-site particles (atoms) with ν0 ) 1. Thus, for any pair combination (0β) involving the wall species, the RISM equation eq 2.1 can be strongly simplified. It follows from eq 2.4a and eq 2.4c that ω ˆ (0β) ) 0

β ) 1..K

ω ˆ (00) ) 1

(2.5a) (2.5b)

13388 J. Phys. Chem. B, Vol. 111, No. 47, 2007

Woelki et al. letting the vector b′ r pass through the total volume. In Figure 1b, the radius of the wall particle 0 has grown to infinity so that its surface forms a plane solid/liquid interface that has been located in the x-y plane of a Cartesian coordinate system. The negative half space -∞ < z < 0 corresponds to the solid phase while the liquid spans the positive half space 0 < z < +∞. Using cylinder coordinates, the volume integration in eq 2.8 +∞ +∞ r ) 2π∫-∞ ∫|z-t|s ds dt and, with eq can be written as ∫V d3 b′ 2.9a,b, eq 2.8 becomes K

hβ(z) )

K

hˆ (0β) )

∑cˆ

(0κ)


β ) 1..K

(2.6)

κ)0

With eq 2.6, we immediately carry out step 2 by forming the limit F0 f 0. This yields K

hˆ (0β) )

∑cˆ

(0κ)


β ) 1..K

(2.7a)

κ)1

Note that a fully equivalent, but not identical, relation is obtained if the above derivation leading to eq 2.7a is repeated for pair combinations (β0) instead of (0β): K

hˆ (β0) ) hˆ (0β) )

∑ωˆ

(cˆ (0λ) +

λ)1

∑ν F

m (λκ) (0κ) ˆ hˆ ) κ κc

κ)1

Both equations, eq 2.7a and eq 2.7b, can be used equivalently to study phenomena occurring in infinite dilution of a monoatomic component 0 (e.g., hydration of simple ions58). However, the structure of eq 5 is more complicated: It involves triple products ω ˆ (βλ)cˆ (λκ)hˆ (0κ) corresponding to complicated double convolutions in b r space while eq 2.7a implies single convolutions only. Hence, eq 2.7b is a more suitable candidate for the limiting process r0 f ∞ in the remaining step 3. The inverse transform of eq 2.7a is K

h(0β)(r) )

∑ ∫c κ)1

V

(0κ)

(r′)(ω(βκ)(|b r - b′|) r +

(βκ) ν κF m (|b r - b′|)) r d3 b′ r κh

β ) 1..K (2.8)

The geometrical meanings of the vectors b r and b′ r are shown in Figure 1a. In b r space, the intramolecular correlation functions ω(βκ)(r) read ω

(r) ) 0 β and κ do not form part of a common molecule (2.9a)

(βκ)

νκ

ω

(r) )

(βκ)

∑

δ(r - l

κ)1

β*κu

)

4π(l β*κu)2 β and κ do form part of a common molecule (2.9b)

u)1

Figure 1a,b illustrates the limiting process r0 f ∞. In Figure 1a, all fluid interaction sites, including the wall particle 0, have a definite radius, and the volume integration is carried out by

c (t)

νκ

∞

|z - t|

1

∑

u)12(l

β*κu 2

)

h(βκ)(s)s ds dt +

∫

+∞

|z - t|

sδ(s - lβ*κu) ds dt (2.10)

The pair correlation functions involving the wall particle 0 (h(0β)(r) and c(0β)(r) in eq 2.8) are now written as singlet correlation functions hβ(z) and cβ(z). With respect to eq 2.9a, the index β at the outer sum in the second term on the righthand side of eq 2.10 indicates that the sum is taken only over sites κ belonging to the same molecular species as site β. The remaining pair correlations h(βκ) have in no way been affected by the limiting process r0 f ∞ and are still bulk total correlation functions. Thus, the inner integrals of the first term on the right-hand side of eq 2.10 can be calculated from the bulk correlations. (V) ) 2π

(βκ)

H

∫

+∞ (βκ)

h

|V|

(s)s ds

(2.11)

The inner integrals of the second term of eq 2.10 are

∫

β ) 1..K (2.7b)

+∞ κ c (t) -∞

∫

+∞ κ

-∞

∑∫

K

(βλ)

m κ κ

κ)1 K β

Figure 1. Geometrical meaning of integration variables. (a) Homogeneous (bulk) site-site system, eq 6; (b) inhomogeneous wall-fluid site system, eq 2.10.

Inserting eq 2.5a,b into eq 2.1, we obtain

∑2πν F ∫

+∞

sδ(s - l β*κu) ds ) |z - t|

{

0 l β*κu

|z - t| > lβ*κu |z - t| < lβ*κu

(2.12)

It is seen from eq 2.12 that this integral is nonzero only in the range (z - lβ*κu) < t < (z + lβ*κu). Thus, the outer integration +∞ z+lβ/κu in the second term of eq 2.10 can be replaced by ∫z-l ∫-∞ β/κu and eq 2.10 can be rewritten as K

hβ(z) )

∑ν F ∫ m κ κ

+∞ κ

-∞

(z - t) dt +

(βκ)

c (t)H

κ)1

νκ

K

∑∑ β

1

κ)1 u)12l

β*κu

∫

z + lβ*κu κ

z - lβ*κu

c (t) dt

β ) 1..K (2.13)

To avoid numerical singularities, the case β/ ) κu with l β*κu ) 0 in the second term on the right-hand side of eq 2.13 is treated separately by forming the limit for l β*κu f 0. 1 f 02l β*κu

lim

lβ*κu

∫

z + lβ*κu κ

z - lβ*κu

c (t) dt ) cβ(z)

(2.14)

With eq 2.14 and the definition of the indirect singlet correlation function γβ ) hβ - cβ, we finally find K

γβ(z) )

∑ν F ∫

m +∞ κ (βκ) (z κ κ -∞ c (t)H

- t) dt +

κ)1

νκ

K

∑∑ β

1

κ)1 u)1 2l κu*β/

β*κu

∫

z + lβ*κu κ

z - lβ*κu

c (t) dt

β ) 1..K (2.15)

Equation 2.15 is a general singlet RISM integral equation (hereafter referred to as SRISM equation) for fluids composed

Singlet-RISM Theory for Solid/Liquid Interfaces

J. Phys. Chem. B, Vol. 111, No. 47, 2007 13389

of one or more molecular species in the vicinity of a plane impenetrable wall. As input, the complete set of (bulk) sitesite total correlation functions h(Rβ) is required which define the quantities H(Rβ)(r) (see eq 2.11). Equation 2.15 is exact and does not contain any approximations. However, the system of eqs 2.15 can only be solved together with a second set of equations relating γβ and cβ, the closure relations. The latter are always approximate. Closure Relations. For inhomogenous solid/liquid systems, the closure relation can generally be written as cβ(z) ) exp[-uβ(z)/kT + γβ(z) + Bβ(z)] - γβ(z) - 1 z > 0 (2.16) z 0 (2.17b)

so that for γ(z) - uβ(z)/kT > 0 the KH closure simply corresponds to the MS approximation cβKH(z) ) -uβ(z)/kT. Third, the modified Duh-Haymet (MDH) closure79 uses a semiempirical bridge function BMDH known to yield quite accurate results for simple Lennard-Jones (bulk) fluids and their mixtures. It is given by β BMDH (z) )

∑F cˆ κ)1

(2.17a)

Second, we employ the so-called partially linearized hypernetted chain or Kovalenko-Hirata (PLHNC or KH) closure59 which is an optimized coupling of the HNC and the mean spherical (MS) approximation. The KH closure has proven quite successful in the description of pure liquids and mixtures in both homogeneous59-63 and heterogeneous72-75,77,78 systems within the framework of RISM and 3D-RISM theories. The KH bridge function reads

{

K

hˆ (0β) ) cˆ (0β) +

K

where T is the temperature, k is Boltzmann’s constant, and uβ(z) is the interaction potential between the wall and an interaction site of group β. The exact bridge functions Bβ(z) remain intractable and are subject to specific approximations. In this work, we check three approximations for their applicability within the framework of the SRISM theory: The well-known hypernetted chain (HNC) closure approximates the bridge function Bβ(z) simply by

BβKH(z)

metrical forces. For such systems, both eq 2.7a and eq 2.7b simplify to ordinary Ornstein-Zernike (OZ) equations:

{

-0.5(γ β)2(z)/(1 + 0.8γ β(z)) γ(z) g 0 γ(z) < 0 c -0.5(γ β(z))2

(2.17c)

Following the usual naming convention, we refer to the combination of SRISM equations and HNC, KH, or MDH closures as SRISM/HNC SRISM/KH and SRISM/MDH theory, respectively. Simple Liquids. In the previous section, we came across two different but equivalent RISM equations, eq 2.7a and eq 2.7b, for correlations involving an infinitely diluted spherical singlesite solute and a molecular solvent. Because of its more complicated structure, eq 2.7b has been rejected and eq 2.7a has been used for the subsequent derivation of the SRISM equation, eq 2.15. Interestingly, the situation is different for simple liquids composed of spherical particles interacting via central sym-

With eq 2.21, the system of eq 2.19b can be rewritten as K

γ β(z) )

∑F ∫ κ)1

+∞ κ (κβ) (z κ 0 h (t)C

- t) dt + u βc (z)

β ) 1..K (2.22)

Equation 2.22 is usually referred to as the SOZ equation68-70 per se while, for obvious reasons, eq 2.19a has never been employed to inhomogeneous systems. Inside the Wall: the Negative Half Space. We have shown in the previous sections that the solution of the SRISM equation, eq 2.15, involves the calculation of the singlet direct wall-fluid site correlation functions cβ(z) in the negative half space -∞ < z < 0, although no fluid particles are located in this region. However, even in the bulk phase, a clear physical interpretation of direct correlations is not at hand. Therefore, some preliminary considerations about the behavior of direct correlations in the negative half space may be useful. A convergent solution of eq 2.15 requires the existence of limits cβ+ and cβ-, both in the positive and in the negative half space. In analogy to site-site correlations in the bulk phase, the former is simply given by cβ+ ) lim c β(z) ) 0 zf+∞

(2.23)


Woelki et al.

To gain some information about cβ- ) limzf-∞cβ(z), we form the limit z f -∞ of eq 2.15. With the substitution p ) z - t and limzf-∞ γ ) -1 - c-, we obtain K

∑ν F

m κ κ

cκ-

κ)1

∫

+∞

H

-∞

(p) dp + cβ- +

(βκ)

νκ

K

∑ ∑c β

κ -

) -1

the discrete lattice of LJ spheres to a smooth plane with homogeneous surface particle density ξ0. Thus, the sum of individual interactions is approximated by an integral. Obviously, the accuracy of the “smooth wall” approximation increases with increasing ratio of fluid site to wall particle size. The Lennard-Jones interaction u(0β) LJ between a single wall particle 0 and a fluid site β at a distance r is

β ) 1..K (2.24) u(0β) LJ (r) ) 40β

κ)1 u)1 κu*β/

Equation 2.24 can be strongly simplified. With eq 2.11, the integrals in the first term on the left-hand side of eq 2.24 are νκFm κ

∫

+∞

-∞

H

∫ ∫ 2πsh (s) ds dp ) ∫ 4πs h (s) ds ) N (2.25a)

(p) dp ) νκFm κ

(βκ)

+∞

-∞ +∞

νκFm κ 0

+∞

(κβ) ex

where N(κβ) ex is the bulk intermolecular excess of sites β in the vicinity of a site κ or, in other words, the bulk intermolecular excess coordination number of a site κ with respect to sites β. Furthermore, the second and third term on the left-hand side of eq 2.24 can be rearranged to give νκ

K

cβ- +

∑ ∑c β

K

κ -

) cβ- + (νβ - 1)cβ- +

κ)1 u)1 κu*β/

∑νc β

κ)1 κ*β

κ κ -

)

) -1

β ) 1..K

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

(2.26)

3. Wall-Fluid Site Interaction Potentials In nonprimitive computer simulations of inhomogeneous fluids, the wall is frequently modeled as a monolayer of fixedin-space Lennard-Jones spheres.15-17,19 For such systems, the total interaction potential between a given fluid site and the wall is calculated by summing up the interactions of the fluid site with all individual wall spheres. In contrast, the SRISM approach models the wall as an infinitely large but structureless sphere. For comparison with inhomogenous simulation results, we mimic the structured wall of the simulations by smearing out

6

(3.1)

(3.2)

∫u z

∞ (0β) LJ (r)r

dr

(3.3)

[ ( ) ( )] 1 σ0β 5 z

10

-

1 σ0β 2 z

4

(3.4)

However, in many systems of practical interest, the wall is not a monolayer but rather an extended solid phase. For such systems, the infinite plane can be replaced by a half space of particle density F0 extending from -∞ < z < 0. In this case, the total wall-fluid site interaction potential is given by

∫

uβp (z) ) 2πF0

κ)1

As before, the subscript β at the square bracket indicates that the term in brackets is added only if sites κ and β form part of a common molecule including the case κ ) β. Remember that νκ is the number of equivalent sites of species κ within a given molecule. Thus, the term [+νκ]β in eq 2.26 is the intramolecular contribution to the total excess coordination number N (βκ) tot . For reasons of material balance, all sites β of a given nonatomic molecule are equally coordinated by a site of group κ (i.e., have identical values of N (βκ) tot ). Thus, all equations of the linear system eq 2.26 with representative sites β being part of such a molecule are identical so that the coefficient matrix of eq 2.26 is singular. Therefore, in our case of aqueous solutions, eq 2.26 does not allow us to calculate the individual limits cβfor the water interaction sites but rather yields a sum rule for the water molecule that must be satisfied by any consistent solution of the SRISM eq 2.15. The singularity of eq 2.15 is of consequence for the convergence of numerical solutions of eq 2.26. We will return to this point later.

σ0β r

With eq 3.1, the integration yields a 10-4 potential

κ)1

K

(βκ) (βκ) κ cex [+νκ]β)Ntot

uβm(z) ) 2πξ0

κ κ -

(2.25b)

-

The total interaction potential uβm(z) between a plane monolayer of wall particles and a fluid site of group β at a distance z is calculated from

∑νc β

12

1 σ0β ) (σ00 + σββ) 2

0β ) (00‚ββ)1/2

K

Inserting eqs 2.25a,b into eq 2.24, we find the limit z f -∞ of the SRISM equation eq 2.15 simplifies to

∑(N

σ0β r

The mixed Lennard-Jones parameters 0β and σ0β are calculated from the Lorenz-Berthelot mixing rules:

(βκ)

|p| 2 (βκ)

[( ) ( ) ]

z

+∞ (0β) uLJ (r)(r

- z)r dr

(3.5)

With eq 3.1, we obtain the 9-3 potential

[ ( ) ( )]

uβp (z) ) 4πF00Rσ30R

1 σ0R 45 z

9

-

1 σ0R 6 z

3

(3.6)

The 9-3 potential uβp (z) is given for completeness only. For comparison with computer simulation results, we have exclusively employed the 10-4 potential, eq 3.1. 4. Models We have solved the SRISM equation, eq 2.15, in conjunction with the HNC closure, eq 2.17a, the KH closure, eq 2.17b, and the MDH closure, eq 2.17c, for pure water and two aqueous electrolyte solutions with salt concentrations 0.25 M and 1 M in the neighborhood of a plane impenetrable wall. For comparison with existing inhomogenous computer simulation data, the system parameters are adapted as close as possible to the inhomogenous systems investigated by Crozier et al.17 using MD simulation methods. Unfortunately, because of the complete lack of repulsive forces acting on the SPC/E water80 hydrogen sites, the RISM approach cannot deal with the specific water model employed by Crozier et al.17 Thus, we use a modified, RISM-adapted SPC/E water model, proposed by Yoshida et al.62 and denoted as SPC/Ef in the following, which differs from the original SPC/E model in the assignment of nonzero LennardJones parameters σHH ) 1.0 Å and HH ) 3.196 × 10-22 J to the hydrogen sites. The ions and the wall atoms are modeled analogously to the work of Crozier et al.17 The interaction parameters and number densities used in the calculations are



TABLE 1: Number of Equivalent Atoms ν, Lennard-Jones Parameters σ and E, Charges q, and Number Densities Gm for the Interaction Site Groups of the Investigated Aqueous Electrolyte Solutions with Salt Concentrations 0.00 M, 0.25 M, and 1.00 M interaction site groups oxygen

hydrogen SPC/E

ν σ/Å /10-22 J q/eo Fm/nm-3 (0.00 M) Fm/nm-3 (0.25 M) Fm/nm-3 (1.00 M) ξ0/nm-2 a

1 3.166b 10.794b -0.8476b 33.3400 33.3400 33.3400

2 0.000b 0.000b

+0.4238b 33.3400 33.3400 33.3400

1.000c 3.196c

anion

1 3.166d 10.793d +1.0000d 0.0000 0.1505 0.6022

1 3.166d 10.793d -1.0000d 0.0000 0.1505 0.6022

wall atom

1.500d 6.903d 0.0000d

44.4444d

a ξ is the surface number density of the wall atoms. b Original Values for SPC/E Water.80 0 from Crozier et al.17

summarized in Table 1. Note that the ionic parameters for the cations and anions are identical except for the electrical charges. 5. Computational Details Bulk Site-Site Correlations. As shown in the previous sections, the SRISM calculation of the singlet density distributions between wall and fluid sites requires the input of corresponding bulk site-site total correlation functions. The latter are calculated from the DRISM/KH theory51,53 on an equidistant grid of 213 sampling points with distance 2 × 10-2 Å. The influence of salt concentration on the dielectric constant of the solution is neglected, and we use a target relative permittivity of rel ) 78.54 at T ) 298.15 K throughout. The system of nonlinear integral equations is solved by an extension of the method by Labik et al.81 Additionally, to estimate the effects of SPC/E versus SPC/ Ef water models and the accuracy of the DRISM/KH bulk correlation functions in general, we have carried out two Monte Carlo computer simulations for the 1 M bulk solution: one with the original SPC/E water model and a second one with the SPC/ Ef water model as the solvent. The 1 M electrolyte solution was formed from 886 water molecules and 16 cations and anions in a cubic simulation box with edge length 29.84 Å. We have used the standard NVT method with Metropolis sampling and periodic boundary conditions, where the intermolecular shortrange interactions are spherically truncated at half of the box length. Long-range Coulombic interactions are treated with the Ewald summation method. After equilibration, each MC run was perfomed over approximately 18 × 106 configurations. From the ion-oxygen pair correlation functions g(+O)(z) and g(-O)(z), the ionic hydration numbers N(O can be calculated: N(O ) 4πFm O

cation SPC/Ef

∫

d(O min ((O)

r)0

g

(r)r2 dr

(5.1)

where d(O min is the position of the first distinct minimum in g((O)(r). Singlet Wall-fluid Site Correlation Functions. Using the DRISM/KH bulk total correlation functions h(βκ)(r), we first calculate the functions H(βκ)(V) from eq 2.11 by numerical integration. Then, for the calculation of the direct and indirect singlet correlation functions cβ(z) and γβ(z), the SRISM/HNC, SRISM/KH, and SRISM/MDH integral equations, eq 2.15 with eq 2.16 and eqs 2.17a-c, are discretized on a non-equidistant grid of N ) 799 sampling points extending from the boundary zmin (in the negative half space) to the boundary zmax (in the positive half space). The grid distance is proportional to |z3| and the value of zmax is fixed to + 81.90 Å. As shown in section 2, the system of SRISM equations becomes more and more ill-

c

Values for RISM Adapted SPC/Ef Water.62

d

Values

conditioned with increasing negative distance from the wall. Thus, the value of the negative grid limit zmin is subject to certain restrictions: On the one hand, it must be large enough to properly include the run of the direct correlation functions toward their negative limiting values; on the other hand, it must be small enough to avoid numerical instabilities due to poor conditioning at large negative distances (see eq 2.26 and subsequent remarks). We have found that, for the systems under consideration, a value of zmin ) -13.06 Å is a good tradeoff. This point will be discussed in more detail in the following section. The grid distance has a minimum of 0.02 Å at z ) 0 and symmetrically increases in both directions. The outermost step widths are 0.14 Å in the negative and 0.60 Å in the positive half space. The resulting system of nonlinear equations is solved iteratively by means of the Newton-Raphson method82 where the Jacobian is obtained as algebraic expression from eq 2.15 and eq 2.16 with eq 2.17a, eq 2.17b, or eq 2.17c, respectively. For the systems under consideration, a root-mean-square residual R e 1 × 10-6 was reached within 4 to 5 Newton steps from initial estimates γβ(z) ) 0 for all components β ) 1..K. Every Newton step requires inversion of the NK × NK Jacobian (1598 × 1598 for pure water and 3196 × 3196 for the electrolyte solutions) which, with direct methods such as LU-decomposition, is quite time consuming. To accelerate this step, we have employed a very fast and efficient iterative algorithm where the inversion of the Jacobian is carried out by means of a Gauss-Seidel iteration83 with the intermediate solution vector being improved after every 15 iterations by a DIIS step.84 Using this method, the numerical solution of the SRISM/HNC, SRISM/ KH, and SRISM/MDH equations for the systems at hand requires no more than 5 minutes on a usual PC with a 32 bit processor of clock rate 1 GHz and 256 MB RAM. Characteristic Quantities for Solid/Liquid Interfaces. From the indirect and direct singlet correlation functions, γβ(z) and cβ(z), we first calculate the total singlet correlation functions hβ(z) ) γβ(z) + cβ(z) and the wall-fluid site density distributions gβ(z) ) hβ(z) + 1. Furthermore, hβ(z) is used for the calculation of some additional characteristic quantities from the following formulas: surface excess Γβ: Γ β ) Fm β

∫ h (z) dz ∞ β

0

(5.2)

electric charge density Fel(z): K

Fel(z) ) eo

∑z ν F

m β β β β h (z)

β)1

(5.3)


Woelki et al.

electric field vector E(z): Eel(z) ) -

eo o

K

∑z ν F ∫ h (t) dt

β)1

m ∞ β β β β z

(5.4)

electrostatic potential φ(z): φ(z) )

eo o

K

∑z ν F ∫ h (t)(z - t) dt

β)1

m ∞ β β β β z

(5.5)

The numerical integrations in eqs 5.2, 5.4, and 5.5 are carried out by means of the trapezoidal rule. 6. Results and Discussion In this section, we present the results of our SRISM calculations for pure water and two aqueous electrolyte solutions with different salt concentrations adjacent to an impenetrable interface. For comparison, Figures 3-5 and 7-9 also include corresponding results of MD computer simulations by Crozier et al.17 for similiar but not identical inhomogenous systems. As mentioned in the previous sections, the differences particularly pertain to the water models and to the wall-fluid site interaction potentials employed in the simulations and in the SRISM calculations: The “reference” computer simulations17 employ the original SPC/E80 water while we have carried out the SRISM calculations using a modified SPC/Ef62 (see section 4) water model as the solvent. Furthermore, the wall is modeled as a fixed-in-space grid of discrete Lennard-Jones spheres in the reference simulations17 while the SRISM approach deals with a structureless, smooth interface. Although we use a wall-fluid site interaction model that matches the structured wall situation investigated in the reference simulations17 as close as possible (see section 3), the wall-fluid site interaction potentials employed in the simulations and in the SRISM calculations certainly differ at small distances. Obviously, such physical differences do affect both the SRISM and the simulation results rendering a quantitative comparison of the wall-fluid site density distributions somewhat difficult. However, as will be seen in the following sections, these discrepancies show up in a consistent and comprehensible manner. Bulk Correlation Functions. As input, the SRISM approach requires the total radial site-site correlation functions h(Rβ) ) g(Rβ) - 1 of the corresponding homogeneous bulk system. These are shown in Figure 2 for the 1 M electrolyte solution as obtained from the DRISM/KH theory for the SPC/Ef water model (solid lines) and from MC computer simulations for both the SPC/Ef (dashed lines) and the SPC/E (dotted lines) water model. For the salt concentrations 0.25 M and 0.00 M, the solvent-solvent (and ion-solvent) correlations are nearly identical with those of Figure 2, and the ion-ion correlations are quite similiar. For the SPC/Ef solvent, the g(OO), g(OH), and g(HH) DRISM/ KH site-site correlations (solid lines) are in good agreement with the MC simulation results (dashed lines). There are only minor deviations concerning the heights and positions of the peaks. These correlation functions, on the other hand, are somewhat different from those obtained from MC calculations with the original SPC/E model (dotted curve). Interestingly, the additional weak oxygen-hydrogen and hydrogen-hydrogen LJ interactions in the SPC/Ef model (not present in the original SPC/E model) primarily affect the oxygen-oxygen correlations: Compared to the SPC/E water, the first minimum and

Figure 2. Site-site pair correlation functions g(r) for a 1 M aqueous (bulk) electrolyte solution. Full lines are calculated from the DRISM/ KH theory for the RISM-adapted SPC/Ef model62 as solvent. Dashed and dotted lines correspond to the results of Monte Carlo simulations for the SPC/Ef62 and for the original SPC/E80 water model, respectively.

the second maximum of g(OO) are more pronounced and clearly shifted to the outside. The ion-oxygen and ion-hydrogen correlations, g((O) and g((H), depict the hydration of the ions. For both cations and anions, the correlation functions form damped alternating oxygen and hydrogen peaks corresponding to hydration layers around the ions. The DRISM/KH cation-solvent correlations g(+O) and g(+H) are in good agreement with both the SPC/E and the SPC/Ef MC results although the initial peaks are somewhat too low. In contrast, the DRISM/KH anion-solvent correlations g(-O) and g(-H) show severe deviations from both the SPC/E and the SPC/Ef MC simulation data. Moreover, the initial peaks of the DRISM/KH g(-H) and g(-O) correlations imply an intramolecular distance lOH g 1.1 Å which is larger than the fixed input value of lOH ) 1.0 Å. This is a typical unphysical feature of the RISM approach85 and is not observed in the computer simulation results. From the g(+O) and g(-O) correlations, we have calculated the coordination numbers N(O of water molecules around the ions within the first hydration shell from eq 5.1. The values are given in Table 2. While the DRISM/KH cation-water coordination is somewhat too low, the anionwater coordination is much too high and reflects the unsatisfying results of the DRISM/KH theory with regard to anion-solvent correlations. The ergodicity problems seen in the MC ion-ion correlations g(++) and g(- -) for both the SPC/E and the SPC/Ef electrolyte solutions render a quantitative comparison with the corresponding DRISM/KH results difficult. Nevertheless, the simulations show the characteristic depletion of likewise charged ions both for the anions and for the cations. While the DRISM/KH theory reproduces this effect in g(++), it yields a distinct peak in g(- -) at the anion-anion contact distance which implies the formation of anion-anion pairs. Although the formation of pairs of likewise charged ions has also been observed in other DRISM and 3DRISM calculations86,87 as well as in MD simulations,88-91



TABLE 2: Ion-Solvent Coordination Numbers N(O Inside the First Ionic Hydration Shells for the SPC/E80 and SPC/Ef62 Water Model as Obtained from MC Computer Simulations and DRISM/KH Calculationsa water model method d+O min/Å N+O d-O min/Å N-O

SPC/E MC 3.66

SPC/Ef MC 3.66

7.37 3.38

7.27 3.46

6.56

7.30

SPC/Ef DRISM/KH 3.55 6.00 4.05 9.30 r( hyd.

a Electrolyte concentration is 1.00 M. The radii of the hydration shells correspond to the first distinct minima of the corresponding ion-oxygen correlation functions g((O).

we tend to believe that this effect is an artifact of the DRISM approach for the system under consideration. The DRISM/KH cation-anion correlation function g(+ -) again is in reasonable agreement with both simulation results although the initial peak is somewhat too low. The accuracy of the DRISM/KH bulk site-site correlation functions can be summarized as follows: except for the anionsolvent and anion-anion pair correlations, the results are in acceptable agreement with computer simulations for the SPC/E and SPC/Ef electrolyte solutions. We have also checked another RISM-adapted SPC/E water model (with different LJ parameters for the water hydrogen site) proposed by Kovalenko et al.86 using both the DRISM/KH and the DRISM/HNC theory. It turned out, however, that the DRISM/KH theory in combination with the SPC/Ef water model yields the best overall agreement with the computer simulation results. Wall-Fluid Site Density Distributions. Figure 3 compares the density distributions of oxygen sites (upper panel) and hydrogen sites (lower panel) of pure SPC/Ef water as obtained from the SRISM/HNC, SRISM/KH, and SRISM/MDH theories with corresponding MD-simulation results17 for the original SPC/E water. Concerning the oxygen density profile, the simulation shows two distinct peaks at 2.3 Å and 5.1 Å with amplitudes of 4.0 and 1.7, respectively. The SRISM/HNC, SRISM/KH, and SRISM/MDH theories predict a small shift of the first two peaks of about 0.12 Å toward the interface. The HNC and MDH theories overestimate the height of the first peak while the KH theory underestimates the amplitudes of the first and the second peak. Altogether, the SRISM oxygen density profiles are in good agreement with the simulation results, in spite of the underlying differences in surface roughness. In fact, it has been shown by Yeh and Berkowitz13 that the atomistic surface structure of a solid strongly affects the lateral arrangement of water molecules in the adjacent water layer but leaves the average water density distribution perpendicular to the surface practically unaffected. The agreement of the SRISM hydrogen density profiles with the simulation data is not as good as in the case of oxygen. The simulation short-range structure is characterized by two overlapping peaks at 2.3 and 3.1 Å, a small local maximum at 4.1 Å, and a peak at 5.4 Å. Except for the small maximum at 4.1 Å, all peaks are reproduced in the SRISM results, but the widths and amplitudes of the first two peaks are distinctly different from those of the simulation. The SRISM profiles show the typical ice-like structure of few water layers adsorbed to the interface found in X-ray scattering and infrared spectroscopy experiments.92-95 We speculate that the short-range deviations in the SRISM and simulation hydrogen density profiles are mainly due to different wall-hydrogen interaction potentials arising from the underlying water models: The nonzero LJ

Figure 3. Normalized wall-fluid site density profiles gO(z) (water oxygen) and gH(z) (water hydrogen) for pure water. The circles are the MD simulation results17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/MDH, SRISM/HNC, and SRISM/KH results, respectively.

parameters of the SPC/Ef water hydrogen sites (see Table 1) give rise to a short-ranged wall-hydrogen 10-4 dispersion potential of the type of eq 28 while the SPC/E water hydrogen sites (with HH ) OH ) 0) do not interact with the (uncharged) wall particles at all.17 To investigate the effect of the wall-fluid site potential on the interfacial water stucture, we have carried out additional calculations with the 10-4 potential, eq 3.4, and the 9-3 potential, eq 3.6, using varying Lennard-Jones parameters for the wall atoms. Interestingly, the oxygen and the hydrogen density profiles are quite unsensitive to variations in the attractive part of the potential. Evidently, the interfacial water structure results from steric constraints rather than from the type and strength of short-ranged dispersion interactions with the wall particles. This is an important finding with respect to a comparison of our work with other theoretical results for confined liquid water. Kinoshita and Hirata67 have calculated SPC/E water density profiles in the vicinity of a large macrosphere of 30 times the diameter of the water molecule (≈ 8.4 nm) using a 12-6 Lennard-Jones wall-fluid site interaction potential. The authors find that for macrosphere diameters larger than 10 times the solvent diameter the macrosphere-fluid site density distributions become practically independent of the macrosphere size and conclude that macrospheres of this size can be regarded as models of a flat wall. However, in spite of a deeper potential well, the amplitudes of their macrosphere-fluid site correlation functions are much lower than those in our calculations and in the simulations. For example, the first peak of the macrosphereoxygen site density profile reaches a height of only 1.4 while the simulation reaches a value of about 4 and the SRISM/HNC


Figure 4. Normalized wall-ion density profiles g+ (z) (cation) and g- (z) (anion) for a 0.25 M electrolyte solution. The circles are the MD simulation results,17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/MDH, SRISM/HNC, and SRISM/KH results, respectively.

and SRISM/KH calculations of about 5.0 and 3.3, respectively. We think that the convergence of the macrosphere-fluid site pair correlations to the limiting correlation functions of the flat wall is very slow and that the chosen macrosphere diameters are still too small. A more promising approach is provided by the 3DRISM theory. It has been employed to obtain density profiles of water at structured crystalline59,60 and metal61,62 slabs of finite thickness. Similiar to simulation techniques, the 3DRISM method involves the use of periodic boundary conditions and, therefore, is more appropriate for a fluid in a slit pore than for a fluid in front of an isolated wall. Nevertheless, the oxygen and hydrogen density distributions for the crystalline slabs are quite similiar to our results in Figure 3. Note, however, that compared with the SRISM approach presented in this paper, the computational time and effort for solving the 3DRISM approach is considerably higher. Figures 4 and 5 show the wall-ion density distributions of a 0.25 M (Figure 4) and a 1.00 M (Figure 5) electrolyte solution. In full agreement with simulations, the density profiles of the water sites obtained from the SRISM/HNC, SRISM/KH, and SRISM/MDH theories are scarcely affected by the presence of ions and practically identical to those of Figure 3. In contrast to the oxygen and hydrogen density profiles, the simulation density profiles for the ionic species are not smooth but suffer from ergodicity problems indicating insufficient MD-simulation times. The simulations were run for 0.2 ns17 while comparable simulations of inhomogeneous electrolyte solutions are run for 5-15 ns15,20,21 to obtain smooth ionic density profiles. Despite the statistical uncertainties, the simulations clearly show that both cations and anions are completely displaced from the initial

Woelki et al.

Figure 5. Normalized wall-ion density profiles g+ (z) (cation) and g- (z) (anion) for a 1 M electrolyte solution. The circles are the MD simulation results,17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/MDH, SRISM/HNC, and SRISM/KH results, respectively.

TABLE 3: Ionic Surface Excesses Γ( for the Inhomogenous 0.25 M and 1.00 M Electrolyte Solution Where the Values Have Been Calculated from Eq 5.2 for the SRISM/HNC, SRISM/KH, and SRISM/MDH Results salt concentration

0.25 M

1.00 M

Γ(/nm-2 (KH) Γ(/nm-2 (HNC) Γ(/nm-2 (MDH)

-0.12 -0.13 -0.14

-0.48 -0.52 -0.53

water layer (≈ 3 Å). The SRISM theories, in contrast, exhibit initial peaks at ≈ 2.3 Å for the cation and ≈ 2.1 Å for the anion indicating directly adsorbed ions integrated into the initial water layer. As is seen from Table 3, the total ionic surface excess depends on the closure relation and decreases in the order KH-HNC-MDH. Similar ionic adsorption peaks have very recently been reported by Tanimura and co-workers78 from RISM calculations for an uncharged macrosphere immersed in an electrolyte solution. However, the discrepancy between the SRISM ionic density profiles and the simulation results concerning the existence or nonexistence of direct ionic adsorption peaks may have several reasons: On the one hand, the SRISM direct adsorption peaks may simply be an artifact redolent of the well-known failure of ordinary SOZ theories to reproduce the dewetting of hard walls adjacent to simple Lennard-Jones fluids.8 On the other hand, it is known that the detailed shortrange structure of ionic density profiles in the immediate vicinity of an uncharged wall, as obtained from atomistic computer simulations, strongly depends on the specific wall-ion (and wall-solvent) interactions.15 Thus, the smooth wall model employed in the SRISM approach may indeed allow for a certain amount of ionic adsorption while the granular wall used in the



Figure 6. Wall-fluid site direct correlation functions cO(z) (water oxygen), cH(z) (water hydrogen), c+ (z) (cation), and c- (z) (anion) obtained from the SRISM approach for a 1.00 M electrolyte solution. The solid, dashed, and dashed-dotted lines correspond to the SRISM/ MDH, SRISM/HNC, and SRISM/KH results, respectively.

Figure 7. Charge density distribution Fel(z) (upper panel), cs ) 0.25 M (middle panel), and cs ) 1.00 M (lower panel). The circles are the MD simulation results17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/MDH, SRISM/HNC, and SRISM/KH results, respectively.

reference simulations does not. More careful investigations are required to draw more definite conclusions. Wall-Fluid Site Direct Correlation Functions. Figure 6 shows the SRISM/HNC, SRISM/KH, and SRISM/MDH wallfluid site singlet direct correlation functions cO and cH for the water atoms and c+ and c- for the ionic components of the confined 1 M electrolyte solution. Corresponding results for the 0.25 M solution and for pure water are quite similiar and not shown here. Each of the direct correlation functions shows a characteristic minimum and maximum (or minimum and maximum curvature in the case of cH) located in the intervals -1 Å e z e 0 Å and 1.5 Å e z e 2.3 Å, respectively. They are connected by steep but featureless intermediate parts. In the positive half space, the wall-fluid site direct correlations are quite similiar to the bulk site-site direct correlation functions: after the first maximum, some characteristic structure is seen with one or more small local extrema which, after a few angstroms, passes into the asymptotics cβ(z) ≈ - uβ(z)/kT. Identifying the leftmost values at zmin - 13.06 Å of the direct correlations with the limiting values cβ-, the sum rule of eq 2.26 is very well satisfied for all closures. However, the course of the direct correlation functions in the negative half space and the individual values of the cβ-’s not only depend on the wall-fluid site interaction potentials but also depend on the negative grid limit zmin. This is due to the bad conditioning of the system of equations, eq 2.15, in the negative half space (see eq 2.26 and subsequent remarks). To estimate the effect of variations in zmin on the resulting wallfluid site density distributions, we have repeated the above

calculations with (a) zmin ) -3.86 Å and (b) zmin ) -24.43 Å. In case a, the computational time reduces to about 20%. In the negative half space, there are deviations of up to 25% in the direct correlations cβ(z < 0) (where hβ ) -1). In the positive half space, however, the maximum change in both direct and total correlations is less that 1%. In case b, our convergence criterion (remaining residual R e 10-6) is never reached. Additional tests show that the proper choice of the negative cutoff distance zmin is not really a problem. Reliable total correlation functions are obtained within a wide range of zmin values extending from approximately -18 Å to nearly -3 Å (the lower limit of -18 Å is determined by machine precision). Charge Density Profiles. Figure 7 compares the space charge distributions obtained from eq 5.3 for the SRISM/MDH, SRISM/ HNC, and SRISM/KH calculations and for the simulation results. As the solvent concentration is much larger than the electrolyte concentrations, the space charge distribution is mainly determined by the water atoms while ionic contributions are practically negligible. SRISM calculations and simulations yield similiar peak positions but there are larger deviations in the peak amplitudes, especially in the case of the HNC and MDH closure. Electric Field and Electrostatic Potential Profiles. Figure 8 shows the z component E of the electric field calculated from eq 5.4. The electrostatic potential profiles calculated from eq 5.5 are shown in Figure 9. In both cases, the general structure of the SRISM profiles is quite similiar to the simulation results. The values of the surface potential φ(0) are found to be almost independent of the specific closure relation and are in good agreement with the simulation values.


Figure 8. Electric field profiles E(z) in an aqueous electrolyte solution adjacent to a plane, impenetrable wall for the salt concentrations cs ) 0.00 M (upper panel), cs ) 0.25 M (middle panel), and cs ) 1.00 M (lower panel). The circles are the MD simulation results17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/MDH, SRISM/HNC, and SRISM/KH results, respectively.

7. Conclusions In this paper, we present a novel symmetry-reduced integral equation theory, called the singlet RISM (or SRISM) approach, for the description of inhomogenous interaction site fluids adjacent to a plane impenetrable wall. In combination with various closure approximations, we have solved the SRISM equations numerically for pure water and aqueous electrolyte solutions at an uncharged solid/liquid interface. The results have been compared with MD computer simulations for similiar but not identical systems. It turnes out that the SRISM results are in reasonable agreement with the simulation data for the wall-oxygen as well as for the wall-hydrogen density distributions, while there are pronounced deviations in the wall-cation and wall-anion correlations. In part, such deviations can be traced back to differences in the physical systems employed in the simulations and in the SRISM calculations. As to the wall-hydrogen correlations, the deviations are mainly due to the different water models (SPC/E in the simulations and SPC/Ef in the SRISM calculations). The differences in the ion distributions mainly pertain to the existence (SRISM) or nonexistence (simulations) of direct ionic adsorption peaks. The latter are observed in the SRISM results exclusively and probably are either artifacts or arise from the smooth wall model used in the SRISM approach. Characteristic integral quantities of the inhomogenous solid/ liquid systems, such as the electric field and electrostatic potential profiles, turn out to be little affected by the specific

Woelki et al.

Figure 9. Electrostatic potential profiles φ(z) in an aqueous electrolyte solution adjacent to a plane, impenetrable wall for the salt concentrations cs ) 0.00 M (upper panel), cs ) 0.25 M (middle panel), and cs ) 1.00 M (lower panel). The circles are the MD simulation results,17 and the solid, dashed, and dashed-dotted lines correspond to the SRISM/ MDH, SRISM/HNC, and SRISM/KH results, respectively.

short-range structure of the inhomogenous fluid and are wellreproduced by the SRISM approach. In summary, we think that the symmetry-reduced SRISM integral equation theory provides an effective and rather reliable tool for the investigation of inhomogenous liquids and mixtures. In contrast to other approaches, such as the inhomogenous MOZ or OZ theories employing dipolar hard spheres or central-force models for the solvent molecules, the SRISM approach is not restricted to water and aqueous solutions but can be applied to arbitrary inhomogenous interaction site fluids and mixtures. Future applications will include mixed solvents, such as water/ acetonitrile or water/acetone mixtures, as well as molecular and multivalent ions, such as nitrate and sulfate. A general advantage of the method is the inherent symmetry reduction that provides a basis for a fast and numerically stable algorithm. For the systems investigated in this paper, that is, water and aqueous electrolyte solutions, the computational effort takes only a few minutes on an up-to-date personal computer, much less than required for atomistic computer simulations. However, application of the SRISM approach to charged interfaces involving long-ranged Coulombic wall-fluid site interactions entrains additional problems of convergence and is by no means straightforward. Such systems will be tackled in a subsequent paper. Acknowledgment. The authors wish to express their sincere thanks to Professor Douglas Henderson, Prof. Andriy Kovalenko, Dr. Deszo _Boda, and Professor Georg Schmeer for

Singlet-RISM Theory for Solid/Liquid Interfaces valuable discussions. The simulation data provided by Dr. Paul Crozier, Dr. Richard Rowley, and Professor Henderson and control calculations provided by Professor Georg Schmeer are also gratefully acknowledged. References and Notes (1) Gouy, J. Physics 1910, 9, 457. (2) Chapman, D. L. Philos. Mag. 1913, 25, 475. (3) Kohler, H.-H.; Woelki, S. In Coagulation and Flocculation, 2nd ed.; Stechemesser, H., Dobias, B., Eds.; Surfactant Science Series 126; CRC Press: Boca Raton, FL, 2005; pp 43-70. (4) Woelki, S.; Kohler, H.-H. Chem. Phys. 2000, 261, 411. (5) Woelki, S.; Kohler, H.-H. Chem. Phys. 2000, 261, 421. (6) Woelki, S.; Kohler, H.-H. Chem. Phys. 2004, 306, 209. (7) Blum, L.; AdV. Chem. Phys. 1990, 78, 171. (8) Henderson, D. In: Fundamentals of Inhomogenous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; pp 177-199. (9) Blum, L.; Henderson, D. In: Fundamentals of Inhomogenous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; pp 239276. (10) Carnie, S. L.; Torrie, G. M. AdV. Chem. Phys. 1984, 56, 141. (11) Attard, P. AdV. Chem. Phys. 1996, 92, 1. (12) Philpott, M. R.; Glosli, J. N. J. Electroanal. Chem. 1996, 409, 65. (13) Yeh, I.-C.; Berkowitz, M. L. J. Electroanal. Chem. 1998, 450, 313. (14) Spohr, E. Electrochim. Acta 1999, 44, 1697. (15) Dimitrov, D. I.; Raev, N. D. J. Electroanal. Chem. 2000, 486, 1. (16) Crozier, P. S.; Rowley, R. L.; Spohr, E.; Henderson, D. J. Chem. Phys. 2000, 112, 9253. (17) Crozier, P. S.; Rowley, R. L.; Henderson, D. J. Chem. Phys. 2000, 113, 9202. (18) Spohr, E. Solid State Ionics 2002, 150, 1. (19) Guymon, C. G.; Hunsaker, M. L.; Harb, J. N.; Henderson, D.; Rowley, R. L. J. Chem. Phys. 2003, 118, 10195. (20) Sachs, J. N.; Petrache, H. I.; Zuckerman, D. M.; Woolf, T. B. J. Chem. Phys. 2003, 118, 1957. (21) Qiao, R.; Aluru, N. R. Colloids Surf., A 2005, 267, 103. (22) Pashley, R. M.; McGuiggan, P. M.; Ninham, B. W.; Brady, J.; Evans, D. F. J. Phys. Chem. 1986, 90, 1637. (23) Kjellander, R.; Marcˇelja, S.; Pashley, R. M.; Quirk, J. P. J. Phys. Chem. 1988, 92, 6489. (24) Kjellander, R.; Marcˇelja, S.; Quirk, J. P. J. Colloid Interface Sci. 1988, 126, 194. (25) Kjellander, R.; Marcˇelja, S.; Pashley, R. M.; Quirk, J. P. J. Phys. Chem. 1990, 92, 4300. (26) Ke´kicheff, P.; Marcˇelja, S.; Senden, T. J.; Shubin, V. E. J. Chem. Phys. 1993, 99, 6098. (27) Carnie, S. L.; Chan, D. Y. C. J. Chem. Phys. 1980, 73, 2949. (28) Eggebrecht, J. M.; Isbister, D. J.; Rasaiah, J. C. J. Chem. Phys. 1980, 73, 3980. (29) Henderson, D.; Blum, L. J. Chem. Phys. 1981, 74, 1902. (30) Rasaiah, J. C.; Ibister, D. J.; Stell, G. J. Chem. Phys. 1981, 75, 4707. (31) Dong, W.; Rosinberg, M. L.; Perera, A.; Patey, G. N. J. Chem. Phys. 1988, 89, 4994. (32) Torrie, G. M.; Perera, A.; Patey, G. N. Mol. Phys. 1989, 67, 1337. (33) Berard, D. R.; Patey, G. N. J. Chem. Phys. 1992, 97, 4372. (34) Berard, D. R.; Kinoshita, M.; Ye, X.; Patey, G. N. J. Chem. Phys. 1994, 101, 6271. (35) Berard, D. R.; Kinoshita, M.; Ye, X.; Patey, G. N. J. Chem. Phys. 1995, 102, 1024. (36) Berard, D. R.; Kinoshita, M.; Cann, N. M.; Patey, G. N. J. Chem. Phys. 1997, 107, 4719. (37) Yamamoto, M.; Kinoshita, M. Chem. Phys. Lett. 1997, 274, 513. (38) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1988, 88, 7826. (39) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1988, 89, 3285. (40) Torrie, G. M.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1989, 90, 4513. (41) Patey, G. N.; Torrie, G. M. Chemica Scripta 1989, 29A, 39. (42) Vossen, M.; Forstmann, F. J. Chem. Phys. 1994, 101, 2379. (43) Vossen, M.; Forstmann, F. Mol. Phys. 1995, 86, 1493.

J. Phys. Chem. B, Vol. 111, No. 47, 2007 13397 (44) Kra¨mer, A.; Vossen, M.; Forstmann, F. J. Chem. Phys. 1997, 106, 2792. (45) Vossen, M.; Forstmann, F.; Kra¨mer, A. Solid State Ionics 1997, 94, 1. (46) Lemberg, H. L.; Stillinger, F. H. J. Chem. Phys. 1975, 62, 1677. (47) Chandler, D.; Anderson, H. C. J. Chem. Phys. 1972, 57, 1930. (48) Chandler, D. J. Chem. Phys. 1973, 59, 2742. (49) Chandler, D.; McCoy, J. D.; Singer, S. J. J. Chem. Phys. 1986, 85, 5971. (50) Chandler, D.; McCoy, J. D.; Singer, S. J. J. Chem. Phys. 1986, 85, 5977. (51) Perkyns, J. S.; Pettit, B. M. Chem. Phys. Lett. 1992, 190, 626. (52) Hummer, G.; Soumpasis, D. M. Mol. Phys. 1992, 75, 633. (53) Perkyns, J. S.; Pettit, B. M. J. Chem. Phys. 1992, 97, 7656. (54) Hirata, F. Bull. Chem. Soc. Jpn. 1998, 71, 1483. (55) Chong, S.-H.; Hirata, F. J. Phys. Chem. B 1997, 101, 3209. (56) Kinoshita, M.; Hirata, F. J. Chem. Phys. 1997, 106, 5202. (57) Kinoshita, M.; Okamoto, Y.; Hirata, F. J. Chem. Phys. 1997, 107, 1586. (58) Imai, T.; Nomura, H.; Kinoshita, M.; Hirata, F. J. Phys. Chem. B 2002, 106, 7308. (59) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 1009. (60) Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2001, 349, 496. (61) Kovalenko, A.; Hirata, F. J. Theor. Comput. Chem. 2002, 1, 381. (62) Yoshida, K.; Yamaguchi, T.; Kovalenko, A.; Hirata, F. J. Phys. Chem. B 2002, 106, 5042. (63) Omelian, I.; Kovalenko, A.; Hirata, F. J. Theor. Comp. Chem. 2003, 2, 193. (64) Sullivan, D. E.; Barker, R.; Gray, C. G.; Street, W. B.; Gubbins, K. E. Mol. Phys. 1981, 44, 597. (65) Borstnik, B.; Janezic, D. Chem. Phys. 1989, 130, 195. (66) Borstnik, B.; Janezic, D. J. Mol. Liq. 1991, 48, 293. (67) Kinoshita, M.; Hirata, F. J. Chem. Phys. 1996, 104, 8807. (68) Percus, J. K. J. Statist. Phys. 1976, 15, 423. (69) Perram, J. W.; Smith, E. R. Chem. Phys. Lett. 1976, 39, 328. (70) Henderson, D.; Abraham, F. F.; Barker, J. A. Mol. Phys. 1976, 31, 1291. (71) Akiyama, R.; Hirata, F. J. Chem. Phys. 1998, 108, 4904. (72) Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 1998, 290, 237. (73) Shapovalov, V.; Truong, T. N.; Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2000, 320, 186. (74) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095. (75) Kovalenko, A.; Hirata, F. J. Mol. Liq. 2001, 90, 215. (76) Bertagnolli, H.; Hausleithner, I.; Steinhauser, O. Chem. Phys. Lett. 1985, 116, 465. (77) Tanimura, A.; Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2003, 378, 638. (78) Tanimura, A.; Kovalenko, A.; Hirata, F. Langmuir 2007, 23, 1507. (79) Duh, D. M.; Henderson, D. J. Chem. Phys. 1996, 104, 6742. (80) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (81) Labı´k, S.; Malijevsky´, A.; Vonˇka, P. Molec. Phys. 1985, 56, 709. (82) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C, 2nd ed.; Cambridge University Press: Cambridge, England, 1994; pp 379-383. (83) Jeffreys, H.; Jeffreys, B. S. Methods of Mathematical Physics, 3rd ed.; Cambridge University Press: Cambridge, England, 1988; pp 305306. (84) Pulay, P. Chem. Phys. Lett. 1980, 73, 393. (85) Hirata, F.; Kovalenko, A. Molecular theory of solVatation; Kluwer Academic Publications: Dordrecht, Netherlands, 2003; p 238. (86) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10391. (87) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10403. (88) Dang, L. X.; Pettitt, B. M. J. Chem. Phys. 1987, 86, 6560. (89) Dang, L. X.; Pettitt, B. M. J. Am. Chem. Soc. 1987, 109, 5531. (90) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1990, 94, 4303. (91) Guardia, E.; Rey, R.; Padro, J. A. J. Chem. Phys. 1991, 95, 2823. (92) Toney, M. F.; Howard, J. N.; Richter, J.; Borges, G. L.; Gordon, J: G.; Melroy, O. R.; Wiesler, D. G.; Yee, D.; Sorensen, L. B. Surf. Sci. 1995, 335, 326. (93) Hasegawa, T.; Nishijo, J.; Imae, T.; Huo, Q.; Leblanc, R. M. J. Phys. Chem. B 2001, 105, 12056. (94) Ewing, G. E. J. Phys. Chem. B 2004, 108, 15953. (95) Asay, D. B.; Kim, S. H. J. Phys. Chem. B 2005, 109, 16760.

Liquid Interfaces Part I: Uncharged

Recommend Documents