A Simple Phenomenological Fix for the Dielectric Constant within the

22702

J. Phys. Chem. B 2006, 110, 22702-22711

A Simple Phenomenological Fix for the Dielectric Constant within the Reference Interaction Site Model Approach P. Gonza´ lez-Mozuelos Departamento de Fı´sica, CinVestaV, AVenida IPN 2508, Colonia San Pedro Zacatenco, 07360 Me´ xico, D. F., Mexico ReceiVed: July 19, 2006; In Final Form: September 4, 2006

The effective interactions among ions immersed in water are studied by means of the effective pair potentials (EPPs) [J. Chem. Phys. 2002, 117, 6133] obtained after contracting (integrating out) the degrees of freedom of the solvent molecules. The dressed interaction site theory (DIST) leads to a simple way of adjusting the effective dielectric constant of the model solvent to its experimental value at standard conditions. The molecular structure of the solvent is mirrored in the structure of the short-ranged component of the induced EPPs, with noticeable differences between the cases with trivial (ideal gas) and nontrivial (experimental) values of the dielectric constant. The shape of these EPPs remains almost invariant over the whole range of salt concentrations considered here. The asymptotic behavior of the EPP between two macroions obtained after contracting the supporting electrolyte (water molecules plus small ions) is also briefly discussed.

I. Introduction The properties and behavior of simple and complex electrolyte systems are clearly determined by the characteristics of the underlying solvent. For example, the long-ranged electrostatic interactions among the ionic species are strongly modified by the polar nature of the molecules of the solvent, giving an effective dielectric constant s of the medium, whereas the shortranged solvation forces are regulated by the specific arrangement of these molecules in the vicinity of those ions. The understanding at a fundamental level of these mechanisms has a paramount importance in the elucidation of a great variety of phenomena and processes in chemical, biological, and industrial systems. Because computer simulations of electrolyte solutions which explicitly take into account the molecular structure of the solvent are extremely difficult to implement due to the huge number of solvent molecules that should be included, the most viable approaches come then from the implementation of the liquid theory methodology. A very successful approach is based on the reference interaction site model (RISM) integral equation,1-4 which provides a detailed description of electrolyte solutions. This equation, complemented with an appropriate set of closure relations, allows the theoretical determination of all the twobody correlations in a molecular fluid. A valuable insight on the properties of these multicomponent systems may also be attained from an explicit distinction between observable and unobservable components and the subsequent study of the effective interactions among the particles of the first set.5,6 In certain experimental setups, only the properties of a subset of the species in the system are observed directly, so the presence and characteristics of the remaining components should be inferred from their effect on those detectable. The osmotic pressure, for example, is mostly determined by the effective dielectric constant and solvation forces mentioned above. A general and rigorous theoretical description of these effective interactions within the context of RISM, with particular emphasis on systems with electrostatic interactions, has been presented in a previous work,5 referred here as paper I.

The effective pair potentials (EPPs) among the observable particles are defined by their corresponding two-point correlations.5,6 A full knowledge of the microscopic structure of the fluid is needed for the theoretical determination of these EPPs. The description of this structure in a uniform fluid is provided by the total correlation functions hij(r) ≡ gij(r) - 1, where the gij(r) are the corresponding distribution functions, and the RISM equation relates them to the intermolecular direct correlation functions caij(r). This complete description is then contracted to the subset of observable particles generating a RISM equation that includes only the hij(r) among these particles. The EPPs are then the pair potentials needed to balance this equation with its corresponding closures. One of the main results of the dressed interaction site theory (DIST) presented in paper I is that the EPPs in systems with electrostatic interactions have two clearly distinguishable parts, one related to the short-ranged effective interactions, which includes the solvation effects, the other corresponding to the renormalized electrostatic interaction, in which the effective dielectric constant appears. The simplistic application of the RISM scheme yields wrong (ideal gas) values of s. It has been shown that this problem can be fixed by assuming special forms for the bridge functions baij(r) that appear in the closure relations that complement the RISM equation.3,7 The main goal of the present work is the application of the DIST results for the EPPs among the immersed ions to the analysis of the conditions needed to attain a realistic value for s. From this analysis follows a very simple prescription for the bridge and direct correlation functions that produces nontrivial values of this dielectric constant. This prescription includes a fitting parameter that can be adjusted to produce the correct value of the effective dielectric constant. At this level of description, the water molecules have been removed from the picture, leaving us with a system which only contains free ions interacting through the corresponding EPPs. This model system can then be used, for example, to implement a computer simulation study of ionic solutions which includes solvation effects and thus goes beyond the usual primitive model description of such systems. The scheme discussed here is

10.1021/jp0645869 CCC: $33.50 © 2006 American Chemical Society Published on Web 10/14/2006

Phenomenological Fix for the Dielectric Constant

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22703

closely related to the approach first used by Medina-Noyola and McQuarrie8 (and further extended by other authors9-11) to derive the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential starting from the primitive model by contracting (removing) the small ions while keeping the large colloidal particles, the main difference being that they used the standard Ornstein-Zernike equation instead of the RISM equation. The following section provides a succinct description of the RISM approach and the explicit definition of the effective pair potentials. Section III summarizes the main results from paper I about the DIST analysis of these EPPs. The specific application of these results to the case in which only the solvent has been contracted is presented in section IV. The description of the model aqueous solution used for the explicit calculation of s is given in section V. Section VI presents the results for the EPPs among the simple ions obtained after contracting the water molecules, including the details about the approximations and fitting procedure involved in the determination of a realistic value of s. A general outline of some features of the EPPs among charged macromolecules in solution obtained after contracting the supporting electrolyte (water plus small ions) is presented in section VII. Finally, section VIII provides some concluding remarks. Some details of the derivation of the DIST results for the EPPs are discussed in the appendix. II. Effective Interactions in Multicomponent Systems Let us consider a multicomponent molecular fluid in which each molecule can be modeled as a set of interaction sites. In general, this system has Nm molecular components, and each molecule of type R has νR different kinds of interaction sites. Two or more distinct interaction sites within a molecule may be considered of the same kind if they are equivalent under transformations belonging to the symmetry group of this molecule. Thus, if θR is the number density of molecules of type R in the system, and m(R) i is the number of sites of kind i within a molecule of type R, then Fi ) θRm(R) is the total i number density of this species (of interaction sites) in the system, Nm νR. and the total number of species in the system is M ) ∑R)1 Interaction sites located in different types of molecules necessarily belong to different species. The two-point total correlation function between sites of species i and j, separated by a distance r in a uniform system, can be written down as the sum of two terms

hij(r) ) h0ij(r) + haij(r)

(1)

for i, j ) 1, 2,.., M, where the first term corresponds to the intramolecular correlations and the second term represents the intermolecular correlations. The standard approach for the determination of this functions makes use of the RISM equations, which basically define the direct intermolecular correlation functions caij(r) in terms of the total correlation functions, and which in Fourier space take the form1-4 M M

Fih˜ aij(k)Fj )

ω ˜ ip(k) c˜ apr(k) (ω ˜ rj(k) + Frh˜ arj(k)Fj) ∑ ∑ p)1r)1

(2)

for i, j ) 1, 2,.., M, where ω ˜ ij(k) ≡ Fiδij + Fih˜ 0ij(k)Fj. The Fourier transformed functions are related to the corresponding real space functions by

˜f (k) ≡ 4π

∫0

∞

2

dr r j0(kr) f (r)

(3)

where j0(x) ≡ sin(x)/x is the spherical Bessel function of order zero.

Assuming that all the interactions beyond the chemical intramolecular links are pairwise additive, then the RISM equations are complemented by the closure relations

caij(r) ) -βuij(r) + haij(r) - ln(1 + haij(r)) + baij(r)

(4)

for i, j ) 1,..., M, where uij(r) is the pair potential between an interaction site of species i and one of species j, β ≡ (kBT)-1, kB is Boltzmann’s constant, and T is the temperature of the system. Equation 4 can be interpreted as the definition of the intermolecular bridge functions baij(r). Very little is known about these functions, and usually some approximations must be done at this level. The approximation baij(r) ) 0, for example, corresponds to the well-known hypernetted-chain (HNC) closure for the site-site correlations. In general, a further relationship is needed to connect the intramolecular total correlations with the pair potentials and intermolecular total correlations. Here we will take the simplest approach and assume that all the molecules are rigid, in which case the h0ij(r) are already fixed from the beginning. The present theory, however, can be easily extended to systems with flexible and elastic molecules by lifting this last restriction.12 Equation 2 contains, in principle, information about the correlations involving all the M components of the system. In many circumstances, however, only a subset of the components is directly observed in a given experimental setup, and this partial knowledge of the structure is usually interpreted in terms of the EPPs that account for the measured correlations. From the theoretical side, the determination of these EPPs, mediated by the unobservable components, is performed as follows.5 Let us start by rewriting eq 2 in the matrix form

R H a R ) Ω C a Ω + Ω C a R Ha R

(5)

where [R]ij ) Fiδij, [Ha]ij ) h˜ aij(k), [Ca]ij ) c˜ aij(k), and [Ω]ij ) ω ˜ ij(k). If m < M is the number of observable species in the system, then each one of the matrixes in the previous equation can be partitioned in accordance with this distinction between observable and unobservable species:

(

X X X ) XOO XOB BO BB

)

(6)

where X stands for R, Ha, Ca, and Ω. For all this matrixes XT ) X, where the superscript T means transpose. The subscript O corresponds to the subset of m observable components, whereas the subscript B correspond to the subset of M - m background (unobservable) components. Thus, the dimensions of the corresponding matrixes are m × m, (M - m) × (M m), m × (M - m), and (M - m) × m for those with subscripts OO, BB, OB, and BO, respectively. The density matrix R is diagonal, and therefore ROB ) RTBO ) 0; where 0 is the appropriate null matrix. Let us also assume here that there are no chemical links between any of the components of the subset O with any of the components of the subset B, and therefore, ΩOB ) ΩTBO ) 0. To simplify the notation, we will use from now on: ROO ) RO, RBB ) RB, ΩOO ) ΩO, and ΩBB ) ΩB. Thus, eq 5 can now be rewritten as a set of four equations: a a Rγ ) ΩR CRγ Ωγ + RR HRγ

a ΩR CRη Rη Haηγ Rγ ∑ η)O,B

(7)

for R,γ ) O,B. Using the equation corresponding to the subscripts R ) B and γ ) O, to eliminate HBO from the equation corresponding to the subscripts R ) O and γ ) O, we obtain

22704 J. Phys. Chem. B, Vol. 110, No. 45, 2006

Gonza´lez-Mozuelos

the contracted RISM equation aeff a RO HaOO RO ) ΩO Caeff OO ΩO + ΩO COO RO HOO RO (8)

be assumed here that it takes this form in the asymptotic regime where limrf∞ r υ(r) ) 1. In this case,

υ˜ (k) )

where a a -1 - CaBB)-1 CaBO Caeff OO ) COO + COB (ΩB

(9)

is the m × m matrix of effective RISM (intermolecular) direct correlation functions that reproduce the observable total correlation functions. Finally, the effective pair potentials among the observable components, ueff ij (r), are defined by the equation eff a a eff a caeff ij (r) ) -βuij (r) + hij(r) - ln(1 + hij(r)) + bij(r; {u }) (10)

for i, j ) 1,..., m, where we have indicated here that the intermolecular bridge functions should be evaluated as functionals of the EPPs. Combining this last definition with eq 4 yields the expression a aeff βueff ij (r) ) βuij(r) + cij(r) - cij (r) +

baij(r; {ueff}) - baij(r; {u}) (11) for the EPPs in terms of the bare pair potentials. In this contracted description the unobservable particles are thus removed, leaving us with an effective system where the interactions among the remaining particles are given by the EPPs defined in the previous equations. Again, some approximations should be done at this point in order to get closed expressions of baij(r; {ueff}) and baij(r; {u}) in a a terms of caeff ij (r) and cij(r), respectively, and hij(r). The EPPs defined here only provide the correct two-body distribution functions, and in general, it is not to be expected that higherorder correlations would be accurately accounted by them. On the other hand, in the limit of infinite dilution of observable particles (RO ) 0), the EPPs become identical to the potentials of mean force (PMFs) among these observable particles, which are defined as βwij(r) ≡ -ln(1 + haij(r)). III. Dressed Interaction Site Theory for Effective Pair Potentials The main concern of the present work are those fluid systems in which the electrostatic interactions play a central role in the determination of their microscopic and thermodynamic properties. Thus, the pair potential between a site of species i and a site of species j, located at positions r1 and r2 respectively, can be written down in the form

βuij(r12) ) βusij(r12) +

∫ dr3 ∫ dr4 zi(r13) lbυ(r34) zj(r42) (12)

is the short-ranged part of the where rij ≡ |ri - rj|, interaction, lb ≡ e2/4πr0kBT is Bjerrum’s length, e is the proton charge, 0 is the dielectric constant of the vacuum, r is the relative dielectric constant of the background, and zi(r) is the spherical electric charge distribution (in units of e) centered in an interaction site of species i. The Fourier transform of eq 12 takes the form usij(r)

4π f(k) k2

(14)

where the function f(k) is such that limkf0f(k) ) 1. From eq 3, and using the fact that j0(0) ) 1, it is clear that qi ≡ z˜i(0) is the total charge (in units of e) associated with a site of species i. Furthermore, any such system must fulfill the global electroM neutrality condition ∑i)1 Fiqi ) 0 to be thermodynamically stable. From eqs 4 and 13 follows that the functions c˜ aij(k) are then of the form

˜i(k) lbυ˜ (k) z˜j(k) c˜ aij(k) ) c˜ sa ij (k) - z

(15)

for i, j ) 1, 2,.., M, where the functions c˜ sa ij (k) are the Fourier transformed “short-ranged” (see also next section) components of the intermolecular direct correlation functions. In matrix notation this result takes the form

Ca ) Csa - lbυ˜ (k)ZZT

(16)

where the column vector Z has the components [Z]i ) z˜i(k). This general form (which can also be demonstrated directly from the density functional definition of the RISM direct correlation functions13) and the RISM equation itself are the only elements needed in the derivation of the DIST results for the total correlation functions and effective pair potentials.5,14 The main results from this theory concerning the EPPs are given below, and a brief account of their derivation is given in the appendix. Using the partition ZT ) (ZTA, ZTB), introducing into eq 9 the four submatrixes a sa ) CRγ - lbυ˜ (k)ZRZTγ CRγ

(17)

for R, γ ) O, B, and using the definition -1 WB ≡ (ΩB-1 - Csa BB) -1 ) (IB - ΩBCsa BB) ΩB

(18)

where IB is the (M - m) × (M - m) unit matrix, we get the following general expression for the effective direct correlation functions sa sa sa (r)T Caeff ˜ (r)(k) Z(r) OO ) COO + COBWBCBO - lbυ O ZO

(19)

The renormalized charge distributions are given by sa Z(r) O ≡ ZO + COBWBZB

(20)

and the renormalized electrostatic potential is related to the bare electrostatic potential by the recursive relation

κB2(k) (r) υ˜ (k) 4π

υ˜ (r)(k) ) υ˜ (k) - υ˜ (k)

(21)

with the screening function given by

βu˜ ij(k) ) βu˜ sij(k) + z˜i(k) lbυ˜ (k) z˜j(k)

(13)

Although it is usually assumed that the electrostatic potential is of the form υ(r) ) r-1, for the sake of generality it will only

κB2(k) ≡ 4πlbZTBWBZB

(22)

From the combination of eqs 14 and 21 follows that the



renormalized electrostatic potential takes the form

υ˜ (r)(k) )

4π f (k) k2 + f (k) κB2(k)

takes the form

and therefore the asymptotic behavior of υ(r)(r) will, in general, differ from that of υ(r). Similar results are obtained in the dressed molecule theory (DMT) introduced by Ramirez and Kjellander.15 Both schemes, DIST and DMT, are generalizations to molecular fluids of the dressed ion theory (DIT) first developed by Kjellander and Mitchell.16 From eqs 11 and 19 then follows that the EPPs in charged molecular fluids are rigorously of the general form

˜ seff βu˜ eff ˜(r) ˜ (r)(k) z˜(r) ij (k) ) βu ij (k) + z i (k) lbυ j (k)

(24)

(r) for i, j ) 1, 2,..., m, where z˜(r) i (k) ) [ZO ]i. The renormalized quantities that appear here, however, are not exactly the same that occur in the corresponding DIST expressions for the correlation functions, unless RO ) 0 (see the appendix). It is also important to point out that the separation between shortranged and long-ranged components in the right-hand side of this last equations is not necessarily meaningful in real space. There are some instances in which the determinant

∆B(k) ≡ det(WB-1) ) det(ΩB-1 - Csa BB)

(25)

is equal to zero for certain real finite values of k, whereas the determinant det(ΩB1- - CaBB) is not null for these same values of k. Hence, both terms on the right-hand side of eq 24 have components with purely oscillatory asymptotic behavior at angular frequencies given by those same values of k, but these oscillatory components cancel out each other from the final result for βueff ij (r). Therefore, the asymptotic behavior of this effective potential, which is the main concern of the following sections, will come in general from the remaining part of the electrostatic term in eq 24. IV. Contraction of the Solvent The central interest of this work is the effect of the molecular structure of the solvent on the effective interactions among the ionic components. To address this issue, one should start by contracting only the solvent, that is, by considering that only the components of the solvent molecules belong to the subset B of background particles, and that all the other components (small ions and polyions) belong to the subset O of observable particles. Usually the solvent is composed by just one type of molecule, but in the most general case it may be constituted by a mixture of different chemical species. The only condition is that all the molecules in the solvent are electrically neutral, that is, each molecule has a null electric monopole. In that case, we have

lim ΩBZB ) 0 kf0

βueff ij (r) ≈

(23)

(26)

where 0 is here the appropriate null column vector. Thus, ΩBZB ≈ k2VB for very small values of k (unless all the solvent molecules also have a null electric dipole, in which case we have a higher even power of k), where VB is a (M - m) column vector. From this fact and the results of the previous section (eqs 18 to 24) follows that z˜(r) ˜i(0) ) qi, κ2B(0) ) 0, and i (0) ) z therefore the asymptotic behavior of the EPPs among the ions

where

lb q q s r i j

(

s ≡ lim 1 + kf0

)

κB2(k) k2

(27)

(28)

is the effective (macroscopic) dielectric constant of the solvent. As commented before, unless some special considerations are taken into account, the RISM approach yields wrong values for this effective dielectric constant. The DIST approach discussed in the previous section can be used to clarify the nature of these considerations. Let us start by rewriting the function κ2B(k) defined in eq 22 in the form

κB2(k) ) 4πlbZTBΩBZB + 4πlbZTBΩBCsa BB(IB -1 ΩBCsa BB) ΩBZB (29)

where use has been made of eq 18. Notice that the factor ΩBZB (≈ k2VB for small values of k) appears twice in the second term of the right-hand side, and only once in the first term. Hence, if one just assumes that the functions c˜ sa ij (k) are analytical at k ) 0 (i.e., |c˜ sa ij (0)| < ∞), then the macroscopic dielectric constant s becomes indeed equal to the trivial ideal gas value

s(0) ≡ 1 + 4πlblim

ZTBΩBZB

kf0

k2

4πlb )1+ 3

∑R θRµR2

(30)

where the sum is over all the solvent chemical species, θR is the number density of species R, and µR is the dipole moment (in units of e) of each molecule of that species. In the case of the usual water models, for example, this ideal gas value is about four times smaller than the experimental value w ≈ 80. The EPPs obtained with this ideal gas value were discussed in paper I.5 Nevertheless, eq 29 also provides us with a way around this difficulty. It is possible to get a value of the dielectric constant different from the ideal gas value of eq 30 by positing functions -2 at the origin, that is, that are of the c˜ sa ij (k) that diverge as k form

c˜ sa ij (k) )

4π aij(k) + c˜ mij (k) k2

(31)

where |c˜ mij (0)| < ∞ and |aij(0)| < ∞ (the functions aij(k) are otherwise arbitrary, any change on them compensated by a corresponding change on the truly short-ranged functions c˜ mij (k)). The second term of the right-hand side of eq 29 would then go as k2, and therefore s * s(0), if

lim det ΩB ) 0 kf0

(32)

which is also a sufficient condition for the fulfillment of eq 26. On the other hand, to avoid a mere rescaling of the site charges of the solvent molecules (as proposed by Cummings and Stell17), it is also necessary to require that limkf0det A * 0, where [A]ij ) aij(k).


Gonza´lez-Mozuelos

Looking back at eq 4, and since the total correlations functions should be short-ranged (the system is assumed to be far from the critical point), then it becomes clear that the long-ranged term coming from the c˜ sa ij (k) functions must be canceled out by a similar term coming from the bridge functions, which must then be of the form

b˜ aij(k) )

4π aij(k) + b˜ mij (k) k2

TABLE 1: Lennard-Jones Parameters O H a (Na) b (Cl)

correlation functions take the form

1 j (l k) 2θw 0 HH

h˜ 0HH(k) )

where again < ∞ (and once more, different functions aij(k) correspond to different functions b˜ mij (k)). Thus, eq 4 can now be rewritten as

h˜ 0HO(k) ) ˜h0OH(k) )

(35)

where f (k) has already been introduced in eq 14. Hence, any computer simulation about a given model can, in principle, be fitted with the appropriate choice of the bmij (r). The arguments of this section provide an alternative point of view to the general analysis given by Raineri and Stell about the closure relations that yield nontrivial values for s.7 The emphasis here, however, is on the effective pair potentials instead of the correlation functions, which may allow an extension of previous results to more complex situations. V. Model System Let us consider here a simple model system for an ionic suspension composed by cations, anions (denoted by the subscripts a and b, respectively), and water molecules. These water molecules are in turn constituted by two kinds of interaction sites: one type corresponding to the two hydrogen atoms (denoted by the subscript H), the other one corresponding to the oxygen atom (denoted by the subscript O). A very useful and quite simple model for the water molecules is the SPC/E water model, which has been fitted to reproduce the experimentally measured structure factors of water at standard conditions.19,20 In this model the intramolecular parts of the total

(36)

1 j (l k) θw 0 OH

h˜ 0OO(k) ) 0

(34)

where all the functions are appropriately short-ranged, and therefore |h˜ aij(0)| < ∞. It is important to emphasize here that the previous arguments are exact and general since they follow directly from the rigorous statistical mechanical formulation of the RISM and DIST formalisms. As mentioned before, eq 2 is essentially the definition of the direct correlation functions c˜ aij(k), whereas eq 16 is basically the definition of the functions c˜ sa ij (k). As it is shown in the appendix, the DIST results stated in eqs 18 to 24, and consequently those stated in eqs 28 and 29, are in turn a direct consequence of these two definitions, with no further assumptions or approximations. Hence, the general form of the functions c˜ sa ij (k) given in eq 31 is a necessary condition for s * s(0) in any interaction site model for polar solvents (e. g., those used in computer simulation studies of water.18). Once it has been prescribed a specific form for the, up to now, arbitrary functions aij(k), then eqs 31 and 33 become in turn the definitions of the functions c˜ mij (k) and b˜ mij (k), respectively. To simplify the numerical calculations, it is convenient to assume that these functions are of the form

aij(k) ) f (k) a(0) ij

σi(Å) 3.1656 0.6000 2.3500 4.4000

(33)

|b˜ mij (0)|

haij(r) ) exp(-βusij(r) + haij(r) - cmij (r) + bmij (r)) - 1

βi 0.26066 0.16782 0.21820 0.16780

(37) (38)

where: θw is the number density of water molecules; lHH and lOH are, respectively, the separations between the hydrogen sites and between each hydrogen site and the oxygen site within each water molecule; and use has been made of the relation FH/2 ) FO ) θw. Besides those already indicated in eqs 36 to 38, all the other intramolecular total correlation functions in the simple electrolyte model considered here are identical to zero. The pair potentials among the all interaction sites are of the form

lb βuij(r) ) βuLJ ij (r) + qi qj r

(39)

for i, j ) a, b, O, H, where the site charge distributions have been assumed to be pointlike, and

(( ) ( ) )

βuLJ ij (r) ) 4ij

σij r

12

-

σij r

6

(40)

is the Lennard-Jones pair potential. If this Lennard-Jones potential is taken directly as the short-ranged part of the -1 when r f 0, interaction, then the csa ij (r) will diverge as r complicating the numerical calculations of the correlations in the system. To avoid this difficulty, we will take instead the following partition for the pair potentials:

lb βuij(r) ) βusij(r) + (1 - exp(-Rr))qi qj r

(41)

lb βusij(r) ) βuLJ ij (r) + exp(-Rr)qi qj r

(42)

where

and R is an arbitrary parameter that can be fitted at convenience. As a result,

f (k) ) (1 + k2/R2)-1

(43)

Though this partition of the pair potential introduces a new pole in the renormalized potentials, it does not affect the relevant physical poles closest to the real axis. The specific values of the parameters in the model used here are as follows. Those for the Lennard-Jones parameters are given, at room temperature, by the values in Table 1 and by the Lorentz-Berthelot mixing rules

σij ) (σi + σj)/2

(44)

βij ) (βi βj)1/2

(45)



where H has been adjusted to fit the first peak of the experimental hydrogen-oxygen correlation.20,21 The intramolecular distances within the water molecules are lOH ) 1.0000 Å and lHH ) 1.6331 Å. The charges of the interaction sites are given by qO ) -2qH ) -0.8476, and qa ) -qb ) 1.0000 (thus, Fa ) Fb ) Cs). Bjerrum’s length for vacuum at room temperature is lb ) 557.00 Å (r ) 1). The arbitrary parameter in eqs 41 and 42 is taken as R ) 5 Å-1. Finally, the density of water molecules is θw ) 0.03344 Å-3 at infinite salt dilution (Cs ) 0), and θw ) 0.03272 Å-3 for Cs ) 6.438 × 10-4 Å-3 ) 1.069 M. For intermediate salt concentrations, the density of water molecules is interpolated linearly. These last values, and those corresponding to the Lennard-Jones parameters for the components a and b, have been previously used by Kovalenko and Hirata for modeling Na and Cl ions in aqueous solution.22 VI. Dielectric Constant and Solvation Forces Let us consider now the contraction of the water molecules, that is, species H and O belong to the subset of background species. For convenience, the subscript B of the previous sections will be replaced by w in the present section, the superscript (r) of the renormalized quantities will be replaced by (w), and the (r). The intramolecular resulting EPPs will be denoted by ueff(w) ij structure of the solvent is then given by

Ωw )

(

2 θw(1 + j0(lHHk)) 2 θw j0(lOHk) θw 2 θw j0(lOHk)

)

(46)

and therefore limkf0detΩw ) 0, so limkf0 ΩwZw ) 0 and the inverse matrix in the second term of the right-hand side of eq 29 does not vanish in that limit, even if the c˜ sa ij (k) are of the form given in eq 31. To avoid the situation in which at least part of the correction to the dielectric constant comes from a direct rescaling of the site charges in the water molecules (the prescription previously proposed by Cummings and Stell,17), it is necessary to set to zero at least one of the coefficients a(0) ij , for i, j ) H, O (indeed, a(0) ij is equal to zero if at least one of the subscripts is different from H or O). The simplest approach, though not unique, seems (0) to be just to assume that a(0) HH ) aOO ) 0, so the pathological (0) (0) (0) (0) case aHHaOO ) aHOaOH is avoided, and the number of adjustable parameters is minimized. For convenience, let us assume the form (0) a(0) HO ) aOH ) 4πlbqHqOγ

(47)

for the remaining coefficient. From eqs 28 to 30 follows that under these circumstances the effective dielectric constant of the contracted solvent is given by

s - 1 )

2 ((0) s - 1) 2 - ((0) s - 1)γ

Figure 1. Intermolecular distribution functions among the components of the water molecules in the limit of infinite salt dilution. Comparison of the results corresponding to an ideal gas value of the dielectric constant against those corresponding to a realistic value of the dielectric constant.

(48)

so only γ and (0) s determine the final value of s. Because for the water model used here (0) s ) 19.62, then the experimental value23 s ) 78.21 is attained with γ ) 0.081486. Thus, the corrections to the functions c˜ aHO(k) ) c˜ aOH(k) are comparatively small, and accordingly, as shown below, their effect on the water microstructure is also rather minor. Figure 1 shows the distribution functions gij(r) ≡ 1 + haij(r), for i, j ) H, O, corresponding to the case of γ ) 0 (dotted lines) with those corresponding to the fitting value γ ) 0.081486 (solid lines), at zero salt concentration (i.e., for θw ) 0.03344

Figure 2. The logarithm of the absolute value of the renormalized electrostatic potential obtained after contracting the water molecules times the separation distance, at infinite salt dilution. Comparison of the results corresponding to an ideal gas value of the dielectric constant against those corresponding to a realistic value of the dielectric constant.

Å-3). These results, and all the others reported here, were determined using a HNC-like closure, that is, bmij (r) ) 0, for all the correlations (for γ ) 0 this is strictly the HNC approximation). As it can be seen from this figure, the correction needed to fit the experimental value of the dielectric constant produces only relatively small modifications of the microstructure of the water components, the most noticeable change being the rather minor reduction in the height of the first peaks, although their positions remain basically uncharged. On the other hand, the changes in the effective interactions (and correlations) among any ionic components immersed in this model solvent are radically altered due to the large correction in the effective dielectric constant. This effect is illustrated in Figure 2, which shows the plot of ln(|rυ(w)(r)|) for the cases of γ ) 0 (dotted line) and γ ) 0.081486 (solid line), at θw ) 0.03344 Å-3. In the first case the plot goes asymptotically to -ln(19.62), whereas in the second it goes to -ln(78.21). Thus, the electrostatic interactions at large separation distances are much more weaker in the second case, and therefore much more realistic. It is also clear, on the other hand, that the renormalized electrostatic potential has, in both cases, a strongly oscillatory behavior at shorter distances, with a wavelength of the order of a couple of angstroms, roughly of the size of the water molecules. This short-range behavior, however, is hidden by the short-ranged part of the ueff(w) (r). ij


Gonza´lez-Mozuelos

Figure 3. The short-range components of the effective pair potentials among the immersed ions obtained after contracting the water molecules at infinite salt dilution. Comparison of the results corresponding to an ideal gas value of the dielectric constant against those corresponding to a realistic value of the dielectric constant. Also shown the corresponding Lennard-Jones potentials.

The effects of adjusting the parameter γ are also noticeable precisely in these short-ranged components. To avoid the difficulties described in the discussion around eq 25, for the purpose of the present analysis let us take instead the partition

lb βueff(w) (r) ) βu/ij(r) + q q ij s r i j

(49)

where i, j ) a, b. Thus, the component u/ij(r) includes the solvation effects involved in the effective interactions among the ions in the solution. Figure 3 presents the plots for these components in the cases of γ ) 0 (dotted line) and γ ) 0.081486 (solid line), at zero salt concentration. The plots of the corresponding Lennard-Jones potentials are also shown for comparison. It is clear from this figure that there is an important

effect of the parameter γ on the calculated solvation forces, although in the cases involving two ions with charge of the same sign these effects are masked by the strong electrostatic repulsions between them. In the cation-anion interactions, however, the electrostatic attraction makes the solvation forces much more relevant, and the deep attractive well located at r ≈ 3 Å for the case of γ ) 0 may induce ionic pairing not present in the case with a realistic dielectric constant. At a finite salt concentration of Cs ) 1.069 M, for which θw ) 0.03272 Å-3, the dielectric constant takes the values s ) 19.28 for γ ) 0 and s ) 72.58 for γ ) 0.081486. Both values are smaller than those corresponding to Cs ) 0, as it is to be expected from the reduction in the water density. The change is, however, proportionally more important in the case of γ ) 0.081486 (-7.2%) than in the case of γ ) 0 (-1.7%). The reason for this is that whereas (0) s - 1 is directly proportional to the solvent density, as can be appreciated from eq 30, eq 48 indicates that the dependence of s - 1 on θw is nonlinear. This same equation also shows that, keeping γ fixed, s goes to (0) s when θw goes to zero, which is the correct limit behavior and also justifies the use of this assumption. On the other hand, the short-ranged potentials βu/ij(r) corresponding to this finite salt concentration are practically indistinguishable from the case of infinite salt dilution, so the solvation forces are basically invariant within this concentration range. Thus, the first term on the right-hand side of eq 49 does not change with the increase in salt concentration, whereas the second term changes very little. This stability of the effective potentials ueff(w) (r) can be contrasted to the marked depenij dence of the PMFs with respect to the salt concentration, which is illustrated in Figure 4. Here are shown the cases: Cs ) 0 (solid line), for which s ) 78.21; Cs ) 0.01069 M (dotted line), for which s ) 78.14; and Cs ) 0.1069 M (dashed line), for which s ) 77.61; all three with γ ) 0.081486. The PMFs for the case of infinite salt dilution are identical to the corresponding EPPs given in eq 49. The positions and general shape of the more noticeable features in these PMFs (e. g., the attractive well at r ≈ 2.82 Å in the cation-anion interaction) are then determined by the corresponding features in the EPPs, but for the screening effects due to the finite salt concentrations which induce a faster decay in the PMFs.22,24 These screening effects, characteristic of the asymptotic behavior of the correlations among the immersed ions, are in turn determined by the effective dielectric constant and the solvation forces included in the EPPs. The length scales involved in the screening can also be analyzed through their influence on the effective interactions between two large charged particles immersed in the electrolyte solution. This last approach is discussed in the following section. VII. Electrolyte Mediated Effective Interactions Let us now add a fifth component to the model system considered in the two previous sections. For the sake of simplicity, let us assume that this extra component, macroions denoted by the subscript M, has no intramolecular structure, although this restriction is not essential and can be easily removed to consider more complex situations (as would be the case if one is considering linear polyelectrolyte solutions25). The effective interaction between any pair of these spherical particles of species M is obtained by contracting all the components of the supporting electrolyte, that is, species a and b now also belong to the subset of background species, in addition to the water molecules. Thus, in the present section, the superscript



Figure 5. The logarithm of the absolute value of the renormalized electrostatic potential obtained after contracting the supporting electrolyte times the separation distance, for the case with realistic dielectric constant (γ ) 0.081486), at two different salt concentrations.

Figure 4. Potentials of mean force among the immersed ions for the case with realistic dielectric constant (γ ) 0.081486) at three different salt concentrations.

(se) will be used instead of (r) for the renormalized quantities, the subscript B will be replaced by se, and the effective pair potential will be denoted by ueff(se) MM (r). The asymptotic behavior of this potential is then basically determined by the nontrivial poles of υ˜ (se)(z), that is, by the solutions of the equation

∆se(z)(z2 + f (z)κ2se(z)) ) 0

(50)

closest to the real axis (see eq 23 and the discussion about eq 25). The asymptotic behavior of ueff(se) MM (r) follows from the last term of eq 24, so it takes the form

βueff(se) MM (r) ≈

lb (1)

E r

exp(-η(1)r) a(1)2 M +

lb (2)

E r

exp(-η(2)r) a(2)2 M (51)

for a sufficiently low salt concentration, where iη(1) and iη(2) are the two purely imaginary solutions of eq 50 closest to the

real axis (η(2) g η(1)), the effective dielectric constants E(1) and E(2) come from the corresponding residues of υ˜ (se)(z) (and, in general, differ from the pure solvent value s), and the (l) renormalized charges are given by a(l) ˜(se) M ) z M (iη ). Previous works at the primitive model level (i.e., no explicit solvent molecules) have predicted that at higher salt concentrations these two imaginary solutions merge into a pair of conjugated complex solutions, and consequently the asymptotic behavior of the potential becomes oscillatory.16 For a very low concentration of small ions, however, the form of eq 51 is valid, and in fact η(2) . η(1), so the second term becomes negligible. In this same limit E(1) ≈ s, and the Yukawa form of the well-known DLVO model potential is recovered as the correct asymptotic behavior of the effective pair potential between two macroions. Deviations from this DLVO model may occur at shorter distances due to the short-ranged part of the effective potential and to the overlapping of the renormalized charge distribution,26 but these seem to be significant only in special situations. Since the main concern here are the length scales in the asymptotic form of β ueff(se) MM (r) in the limit of infinite macroion dilution, it is not necessary to specify here the details of the short-ranged interactions involving at least one macroion. Figure 5 illustrates the behavior of ln(|rυ(se)(r)|) for different salt concentrations when FM ) 0. In this limit the asymptotic length scales of υ(se)(r) are the same as those of the correlations and potentials of mean force (illustrated in Figure 4) among the small ions. The plots in Figure 5 correspond to the salt concentrations values Cs ) 0.01069 M ) 6.438 × 10-6 Å-3 (dotted line) and Cs ) 0.1069 M ) 6.438 × 10-5 Å-3 (solid line), again with γ ) 0.081486. The case of Cs ) 0 is identical to the solid line in Figure 2 (both figures are at the same scale to make the comparison easier). As can be seen from these figures, the increase in salt concentration not only changes the slope of the asymptotic part, but also extends the range of the short-ranged oscillatory behavior induced by the molecular microstructure of the solvent. In fact, for Cs ) 1.069 M (not shown here) the Yukawa part is completely hidden by this rapidly decaying oscillatory component, making much more difficult the determination of the salt concentration at which η(2) ) η(1). The determination of the screening parameters η(l) is greatly simplified by going back to the effective pair potentials calculated in the previous section, that is, by considering a simplified model in which the water molecules are not taken explicitly into account and the ions interact through the potential given in eq 49 (i.e., a primitive-like model). This approach takes


Gonza´lez-Mozuelos

advantage of the fact that, for γ ) 0.081486, the potentials u/ij(r) (solid lines in Figure 3) are basically invariant for the whole range of salt concentration considered here, and the change of s with this concentration is easily determined from eq 48 assuming that θw (and consequently s(0) - 1) varies linearly with Cs. For the whole range of salt concentration considered here the value of the screening length η(1) is practically indistinguishable from the standard Debye-Hu¨ckel value κ )

matrix in the last term of this equation can be expanded as

(IB +

lbυ˜ (k)ZBZTBWB)-1

∞

)

(-lbυ˜ (k)ZBZTBWB)n ∑ n)0 ∞

) IB -

lbυ˜ (k)ZBZTBWB

∑×

n)1

x4πlw(Faqa2+Fbqb2) with lw ) 7.122 Å.

(-lbυ˜ (k)ZTBWBZB)n-1 (A2) where the last line follows from taking into account that

VIII. Concluding Remarks

(ZBZTBWB)m+1 ) ZBZTBWB(ZTBWBZB)m The effects of the solvent on the interactions among the immersed ions appear at two different levels. On one hand, there is an effective dielectric constant much more larger than the vacuum value. The primitive model usually applied to the study of electrolyte solutions takes into account this effect by simply assuming such value for the dielectric constant, although its variation with the salt concentration via its effect in θw is typically neglected. On the other hand, the explicit nature of the molecular structure is mostly evident in the short-ranged part of the effective pair potentials. This last effect accounts for the solvation forces among the ions, which are not included in the primitive model. These forces, however, can be easily included in a primitive-like model by replacing the hard sphere potential with the potentials u/ij(r) discussed above. It is possible to attain a nontrivial value of the dielectric constant by sidestepping the assumption that the functions csa ij (r) defined in eq 16 are necessarily short-ranged. This still leaves some leeway in the conceivable forms of these direct correlation functions. The simple approach proposed here reduces this range of possibilities by keeping only the hydrogenoxygen term with an adjustable parameter that is then used to fit the value of s with minimal changes on the correlations among the water components. The advantage of this prescription is that it is easier to implement than the previous approaches suggested by Perkyns and Pettitt3 and Raineri and Stell.7 Furthermore, the EPPs obtained by taking out the solvent molecules can also be used to model the interactions in more complex systems. These EPPs contain all the relevant information about the dielectric constant of the medium and the solvation forces among the ionic species, so they can be used for relatively simple computer simulations and theoretical analyses of the structural and thermodynamic properties of complex electrolyte solutions (that is, without including explicitly the molecular structure of the solvent). For example, the folding of suspended peptides and proteins27-29 and the conformation of DNA in water30-32 can be treated in this way. Also the theoretical study of solutions of charged colloids can be enriched by this approach.26,33-37 Acknowledgment. This work was financially supported by Conacyt-Mexico under Grant SEP-2004-C01-47200/A-1. Appendix A: DIST for EPPs From eqs 17 and 18, we have that

) WB(IB +

Thus,

(ΩB-1 - CaBB)-1 ) WB - lbυ˜ (r)(k)WBZBZTBWB (A4) where υ˜ (r)(k) is the renormalized electrostatic potential defined in eqs 21 and 22. The combination of this last result with eqs 9, 17, and 22 yields the result sa sa sa Caeff OO ) COO + COBWBCBO

( (

-lb υ˜ - υ˜ 2 -lb υ˜ - υ˜

( ) )

κB2 κB2 2 (r) υ˜ ZOZTO + υ˜ 2 4π 4π

)

κB2 (r) sa T υ˜ (ZOZTBWBCsa BO + COBWBZBZO) 4π

T sa -lb υ˜ (r) Csa OBWBZBZBWBCBO

(A5)

from which, after further regrouping of terms and taking into account eq 21, follows eqs 19 and 20. A couple of useful formulas can be obtained from the definition of the functions h˜ sa ij (k) trough the matrix relation

R Hsa R ) Ω Csa Ω + Ω Csa R Hsa R

(A6)

which is analogous to eq 5. From the equivalent of eq 7 then follows that -1 sa sa Csa OBWB ) (ΩO + ROHOORO) ROHOBRB

(A7)

and sa WB ) ΩB + RBHsa BBRB - RBHBORO(ΩO + -1 sa ROHsa OORO) ROHOBRB (A8)

The functions h˜ sa ij (k) are the “short-ranged” components of the corresponding intermolecular total correlation functions, so the two previous equations connect the renormalized quantities υ˜ (r)(k) and z˜(r) i (k) (see eqs 20 and 22) with the global microstructure of the system, and are particularly useful in the limit RO ) 0. References and Notes

(ΩB-1 - CaBB)-1 ) (ΩB-1 - Csa ˜ (k)ZBZTB)-1 BB + lbυ lbυ˜ (k)ZBZTBWB)-1

(A3)

(A1)

where IB is the (M - m) × (M - m) unit matrix. The inverse

(1) (2) 4133. (3) (4) 183.

Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. Hirata, F.; Rossky, P. J.; Pettitt, B. M. J. Chem. Phys. 1983, 78, Perkyns, J.; Pettitt, B. M. J. Chem. Phys. 1992, 97, 7656. Stell, G.; Patey, G. N.; Hoÿe, J. S. AdV. Chem. Phys. 1981, 48,

Phenomenological Fix for the Dielectric Constant (5) Bagatella-Flores, N.; Gonza´lez-Mozuelos, P. J. Chem. Phys. 2002, 117, 6133. (6) Gonza´lez-Mozuelos, P.; Meńdez-Alcaraz, J. M.; Castan˜eda-Priego, R. J. Chem. Phys. 2005, 123, 214907. (7) Raineri, F. O.; Stell, G. J. Phys. Chem. B 2001, 105, 11880. (8) Medina-Noyola, M.; McQuarrie, D. A. J. Chem. Phys. 1980, 73, 6279. (9) Beresford-Smith, B.; Chan, D. Y. C.; Mitchell, D. J. J. Colloid Interface Sci. 1985, 105, 216. (10) Belloni, L. J. Chem. Phys. 1986, 85, 519. (11) Gonza´lez-Mozuelos, P.; Carbajal-Tinoco, M. D. J. Chem. Phys. 1998, 109, 11074. (12) Curro, J. G.; Schweizer, K. S. J. Chem. Phys. 1987, 87, 1842. (13) Sumi, T.; Imai, T.; Hirata, F. J. Chem. Phys. 2001, 115, 6653. (14) Gonza´lez-Mozuelos, P.; Yeom, M. S.; Olvera de la Cruz, M. Eur. Phys. J. E 2005, 16, 167. (15) Ramirez, R.; Kjellander, R. J. Chem. Phys. 2003, 119, 11380. (16) Kjellander, R.; Mitchell, D. J. J. Chem. Phys. 1994, 101, 603. (17) Cummings, P. T.; Stell, G. Mol. Phys. 1981, 44, 529. (18) Saint-Martin, H.; Hernańdez-Cobos, J.; Ortega-Blake, I. J. Chem. Phys. 2005, 122, 224509. (19) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (20) Lombardero, M.; Mart×c1n, C.; Jorge, S.; Lado, F.; Lomba, E. J. Chem. Phys. 1999, 110, 1148. (21) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47.

J. Phys. Chem. B, Vol. 110, No. 45, 2006 22711 (22) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10391; ibid 2000, 10403. (23) Fernańdez, D. P.; Mulev, Y.; Goodwin, A. R. H.; Sengers, J. M. H. L. J. Chem. Phys. Ref. Data 1995, 24, 33. (24) Kinoshita, M.; Harano, Y. Bull. Chem. Soc. Jpn. 2005, 78, 1431. (25) Ermoshkin, A. V.; Kudlay, A. N.; Olvera de la Cruz, M. J. Chem. Phys. 2004, 120, 11930. (26) Carbajal-Tinoco, M. D.; Gonza´lez-Mozuelos, P. J. Chem. Phys. 2002, 117, 2344. (27) Hagihara, Y.; Kataoka, M.; Aimoto, S.; Goto, Y. Biochemistry 1992, 31, 11908. (28) Goto, Y.; Nishikiori, S. J. J. Mol. Biol. 1991, 222, 679. (29) Pliego-Pastrana, P.; Carbajal-Tinoco, M. D. J. Chem. Phys. 2005, 122, 244908. (30) Pohl, F. M.; Jovin, T. M. J. Mol. Biol. 1972, 67, 375. (31) Soumpasis, D. M.; Robert-Nicoud, M.; Jovin, T. M. FEBS Lett. 1987, 213, 341. (32) Huang, C.-I.; Olvera de la Cruz, M. Macromolecules 2002, 35, 976. (33) Hansen, J.-P.; Lo¨wen, H. Annu. ReV. Phys. Chem. 2000, 51, 209. (34) Quesada-Pe´rez, M.; Callejas-Fernańdez, J.; Hidalgo-A Ä lvarez, R. AdV. Colloid Interface Sci. 2002, 95, 295. (35) Aqua, J.-N.; Cornu, F. Phys. ReV. E 2003, 68, 26133. (36) Denton, A. R. Phys. ReV. E 2004, 70, 31404. (37) Anta, J. A. J. Phys.: Condens. Matter 2005, 17, 7935.

A Simple Phenomenological Fix for the Dielectric Constant within the

Recommend Documents