2116
J. Phys. Chem. 1990, 94, 2116-2123
molecules respond to the rapidly alternating electrical field of the vibrating ion pair. This response gives rise to high-frequency librational motions of neighboring solvent molecules. The rotational response is also faster than the translational response. This seems to be in line with a general assumption made in the mean spherical a p p r ~ x i m a t i o n ; ~that ~ . ~is, ~ only the reorientational mechanisms are important in the solvation response. However, the rotational motions of neighboring solvent molecules see a time-dependent inhomogeneous perturbation source. In other words, as far as solvent rotations are concerned, the ion pair is not a sphere. It is difficult to see how any type of continuum theory can correctly portray these fundamental dynamical differences between the various types of solvent motions. The CA method was found to provide a poor representation of the motion of the solvent surrounding the rapidly moving ion pair. This is not surprising, since the CA method was really not developed for handling such problems. Strictly speaking, CA is linked to Brownian motion theory where solvent motion is so rapid (23) Calef, D. F.; Wolynes, P.G. J . Phys. Chem. 1983, 87, 3387. (24) Sumi, H.; Marcus, R. A. J , Chem. Phys. 1986, 84, 4894.
compared with that of the solute that solvent readjustment is effectively instantaneous for every change of solute configuration. It is also valid if the “reaction coordinate” executes only small amplitude excursions from its initial value on the time scale of the process of interest.25 However, the “friction”, which guides the reaction rate, is strongly influenced by solvent dynamics perturbed by ultrafast solute motion^.^ Thus, one must use caution when applying CA to MD calculations of ionic reactions or barrier crossing problems, where the reacting particles are subjected to very strong forces and may suffer large amplitude displacements over very short periods of time. Acknowledgment. Financial support at the PQRL has been shared jointly by the Robert A. Welch Foundation (D-0005, 50% and D-1094, 5 % ) , the National Science Foundation (CHE8611381,39%), and the State of Texas Advanced Research Program (1 306, 6%). Computer time was furnished by the Texas Tech University Academic Computing Center and the Pittsburgh Supercomputing Center. (25) Adelman, S. A. J. Chem. Phys. 1984, 81, 2776.
Single-Ion Activity Coefficients and Structure of Ionic Fluids. Results for the Primitive Model of Electrolyte Solutions Peter Sloth* and Torben S . Smensen Fysisk- Kemisk Institut and Center of Modelling. Nonlinear Dynamics and Irreversible Thermodynamics (MIDIT), Technical University of Denmark, Bldg. 206, DK 2800 Lyngby, Denmark (Received: May 23, 1989)
The hypernetted chain (HNC) approximation has been used to obtain single-ion activity coefficients (SIAC) for a number of primitive model 1:1 and 2:l electrolyte solutions. In this work we investigate the variation in the S A C , as well as in the radial distribution functions, by changing the cation-anion radii ratio, keeping the distance of contact between the cations and anions fixed. The HNC results are compared to available Monte Carlo data. Finally, the ability of the primitive model to predict activity coefficients in aqueous KCI solutions is investigated and the results are compared to those found for a recently proposed extended model.
Introduction
Single-ion activity coefficients (SIAC) are widely used in the theory of electrolyte solutions, but unfortunately such quantities cannot be determined by traditional thermodynamic measurements without the introduction of additional information or assumptions. Although some attempts have been made to evaluate SIAC experimentally by introduction of extrathermodynamic approaches,’T2 it still remains a rather controversial subject in the theory of ionic fluids. From a theoretical statistical mechanical point of view, however, we find that SIAC are well-defined quantities, which-at least in principle-might be determined by computer simulations, if the microscopic interactions between molecules in the systems are k n o ~ n . ~The . ~ intermolecular interactions in real electrolyte solutions are very complex, and strongly simplified models are usually used at present in theoretical investigationsof such systems. The simplest possible physical model that gives a reasonable description of some important features of electrolyte solutions is ( 1 ) Milazzo, G.; Bonciocat, N.; Borda, M . Electrochim. Acta 1975, 21, 349. ( 2 ) Hurlen, T. Electrochim. Acta 1975, 20, 499. (3) (a) Sloth, P.; Ssrensen. T. S. Chem. Phys. Lett. 1988, 143, 140. (b) Svensson, B. R.; Woodward, C. E. Mol. Phys. 1988, 64, 247. (4) (a) Sloth, P.; Ssrensen, T. S. Chem. Phys. Lett. 1988, 146, 452. (b) Ssrensen, T. S.; Jensen, J. 8.; Sloth, P.J . Chem. SOC.,Faraday Trans. 1 1989, 85, 2649.
0022-3654/90/2094-2116$02.50/C
the primitive model (PM), which is defined by the pair potential functions uii(r) =
a,
= qiqj/47tr,
r
< a,
(1)
r L a,
In (1) qi is the charge of species i, t the dielectric permittivity of the pure solvent, and aV the distance of contact between an ion of species i and an ion of species j . We note that one of the most significant shortcomings of the PM is that it treats the solvent as a structureless dielectric continuum. Some more refined models, which treat the solvent on a molecular level, have recently been c o n ~ i d e r e d . ~The . ~ investigation of such models is in an initial stage, and only limited preliminary Monte Carlo (MC) calculations are presently available.6 One should expect, however, the PM to yield correct results for the activity coefficients in the limit of vanishing ionic c~ncentration,~ although the radial distribution functions obtained by the PM might be quite inaccurate at small separations.6 (5) (a) Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1988, 88, 7715. (b) Kusalik, P.G.; Patey, G. N. J. Chem. Phys. 1988,89, 5843. (c) Kusalik, P. G.: Patev. G. N. J. Chem. Phvs. 1988. 89. 1478. (6) Ckan, K. Y.; Gubbins,‘K. E.; Henderson, D.; Blum, L. Mol. Phys. 1989, 66, 299. (7) Kusalik, P.G.; Patey, G. N . J. Chem. Phys. 1987. 86, 5110
0 I990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2117
Single-Ion Activity Coefficients For the time being, the PM is the most well studied model of electrolyte solutions, and Monte Carlo data for single-ion activity coefficients have recently been published for a few PM system^.^ Also, Ebeling and Scherwinski have discussed the application of the PM within the mean spherical approximation (MSA) in order to evaluate SIAC in real electrolyte solutions.8 One of the most successful statistical mechanical approximations so far for systems with long-range forces (e.g., Coulomb forces) is the hypernetted chain (HNC) appro xi ma ti or^.^^^^ Also, the H N C approximation can be improved in a quite straightforward manner.]' A further important feature of the H N C approximation, in this context, is that an explicit, closed equation can be obtained for the SIAC, which is exact within the framework of the H N C t h e ~ r y . ] ~ In , ' ~this paper we present some SIAC and radial distribution functions calculated by the H N C approximation for aqueous 1 :1 and 2: 1 electrolyte solutions by application of the primitive model. The results are compared to available M C data, some of which are freshly taken from our computer. Even if the PM may be insufficient to give correct quantitative results for real electrolyte solutions, we believe that some qualitative insight might be obtained, which is valid also for real systems. In the last section of this paper, we discuss the application of the PM to aqueous KC1 solutions and compare the results to experimental data and to an extended model recently considered by Kusalik and pate^.^ It should be noted that even though this paper is primarily concerned with electrolyte solutions, much of the discussion applies equally well to other ionic fluids, Le., fused salts and classical plasmas. Theory We shall only consider systems where the total potential energy is given as a sum of pair potentials. In the following we furthermore assume-for the sake of simplicity-that the particles posses spherical symmetry. The H N C approximation is then given by the following relation between the radial distribution functions, gij(r), and the direct correlation functions, cij(r) gij(r) = exP[-@uij(r)
+ gij(r)
- 1 - cij(r11
(2)
or cij(r)
= exP[-Ouij(r)
+ yij(r)I
- [1
+ yij(r)l
(2*)
with the definition yij(r) = gij(r) - 1 - cij(r)
In (2), uij(r)is the pair potential function between a particle of species i and a particle of species j and 8 1/ k T , k being the Boltzmann's constant and T the absolute temperature. It is remarked that the function yi,(r) has the useful property of being everywhere continuous, which is not the case for cij(r)and gij(r), for hard-core systems. The functions cij(r)and hij(r) gij(r) 1 are related via the Ornstein-Zernike (OZ) equation Yij(rl2)
E
hij(ri2) - Cij(r12) = CPkJCik(r13) hkj(r32) dr3 k
(3)
where Pk is the particle density of species k and the summation is taken over all species in the system. In the following we sketch the derivation of the basic formula for calculations of SIAC within the framework of the H N C approximation, and its numerical application is discussed in the next section. (8) Ebeling, W.; Scherwinski, K. Z . Phys. Chem. (Leipzig) 1983, 264, 1 . (9) (a) van Leeuwen, J. M. J.; Groeneveld, J.; DeBoer, J. Physica 1959, 25, 792. (b) Meeron, E. J . Math. Phys. 1960, I , 192. ( I O ) HNC calculations for the PM were pioneered by Carley [Carley, D. D. J . Chem. Phys. 1967,46, 37831 and Rasaiah and Friedman [Rasaiah, J. C.; Friedman, H. L. J . Chem. Phys. 1968,48,2742] and have since been used in many theoretical studies on ionic systems. (11) Baquet, R.; Rossky, P. J. J . Chem. Phys. 1983, 79, 1419. (12) Verlet, L.; Levesque, D. Physica 1962, 28, 1124. (13) Hansen, J. P.; Torrie, G. M.; Vieillefosse, P. Phys. Reu. A 1977, 16, 2153.
In passing we note that thermodynamic functions calculated for the primitive model are given in the McMillan-Mayer set of independent ~ariab1es.I~This implies that the chemical potential of the solvent is held fixed rather than the pressure (which is the normal laboratory condition). The activity coefficient, yi, of species i can be given by the eq~ation'~ In yi = ~k P k x ' ~ @ u i k ( r 1gik(rl2,h) 2) dr2 dh
(4)
where the pair potential is multiplied by X whenever particle 1 of species i appears. This particle might be thought of as a "test particle" that is gradually introduced into the system as the parameter X is varied from X = 0 to X = l . The H N C approximation to gjk(r&t) is given by the equation (cf. (2)) gik(rl2,X) = exp[-APUik(r12) + Yik(rlZ~X)l (5)
Substitution of (6) into (4) yields
Integration over h in (7) can be performed explicitly, by application of the OZ e q u a t i ~ n . ' ~The * ' ~ final result is
In Yi = -xPkJCik(r) k
dr + ( l /2)ZPk1hik(r) k Yik(r) dr
(8)
It should be mentioned that mean ionic activity coefficients might be estimated for primitive model electrolytes by application of the so-called zero separation theorem within the framework of the H N C approximation.16 This approach, however, seems only to give reliable results for quite dilute systems. It is interesting to note that (8) yields an expression for In yi that does not depend explicitly on the pair-potential functions, but is a function of structural properties only. Thus if the radial distribution functions are known for a given system by some means one might directly estimate approximate activity coefficients by application of this equation-without explicit detailed knowledge about the forces in the system. This would be an interesting approach, e&, for fused salts, for which rather detailed information about the partial structure factors for real systems is known from neutron diffraction measurements.]' Method of Solution In this section, the numerical solution of the H N C approximation is discussed for the primitive model, involving two different species: one type of cations and one type of anions. The procedure is readily generalized to many-component systems. For such a two-component system, the OZ equation is equivalent to (9a-c) in (3D)Fourier space ( k = lkl)
Tii(k)= -fii(k) + [Fii(k)(l-pjZjj(k))+ ~ ~ ? ~ ~ ~ ( k ) ] (9a) /D(k) (9b) Yij(k ) = -Fij( k ) + Fij( k ) / D(k ) with D ( k ) = [ 1 - p,Zii(k)][ 1 - pjZjj(k)]- ~ ~ p , C ~ ~ ~and ( k i) ,# j (9c) where and f i j ( k )are the Fourier transforms of y i j ( r )and cij(r), (14) (15) (16) (17)
McMillan Jr., W. G.;Mayer, J. E. J . Chem. Phys. 1945, 13, 276. Hill, T. L. Statistical Mechanics; McGraw-Hill: New York, 1956. Sloth, P. J . Chem. Phys. 1988, 89, 5078. Enderby, J. E.; Neilson, G. W. Adu. Phys. 1980, 29, 323.
2118 The Journal of Physical Chemistry, Vol. 94, No. 5, 1990
respectively. The H N C equation can then be solved by the following procedure: I . Make an initial guess for Eij(k). 2 . Calculate T j j ( k )by (9a-c). 3 . Take the inverse transform of y i j ( k )to find y i j ( r ) . 4. Make a new guess for cij(r) by application of ( 2 * ) . 5. Fourier transform cij(r) to give t i j ( k )and check if convergence is achieved. If this is not the case, proceed from step 2 . The long-range Coulomb forces have to be taken carefully into account, if reliable results are to be found for the primitive model. m,la we might define a short-ranged Since cij(r) -@uii(r) as r function ~ * ~ ~by( r )
-
-
c*,,(r) = c J r )
+ /?qlq,/4mr
(for all
r)
(10)
Sloth and Srarensen TABLE I: Values of the Excess Internal Energies, Activity Coefficients, and Osmotic Coefficients Obtained by the HNC Approximation for A Q U ~ ~1:l U SElectrolytes ( B = 1.681) d+ld- - E e x / N k T @ In y+ In y0.0004212 0.0047975
2 3 1 2 3 1 2 3
0.046 24 0.090 95
~ * ~ ~= yij(r) ( r ) - Pqiqj/4atr (for all r )
-
m,
(12)
In the HNC approximation - ~ * ~ ~ ( r )is/ fjust l the potential of mean force between an ion of species i and an ion of species j for r 2 aij. The function yij(r)might then be calculated from Tij(k) by the formula yi,(r) = ( 1 / 2 x z r )
1-
T * j j ( k ) ksin ( k r ) dk
+ Pqiq,/4atr
(13)
where T*l,(k) = T J k ) - Pq,q,/ekZ
(14)
In our calculations, the initial guess for tlJ(k)was simply made by setting c*,,(r) = 0 in (1 I ) . The integrations in ( 1 1 ) and ( 1 3 ) were done numerically by application of Filon's method,I9 which has been found to be very satisfactory for trigonometric integrals.20 It is also noted that, for the parameters used in the present calculations, no divergence problems were found by the method outlined above. Thus it was not necessary to mix successive guesses in the iterative procedure, although this might be the case at higher densities. For each system, the integration parameters were adjusted in order to give the optimal numerical performance, and as many as 8001 grid points were used to represent the oscillating t J k ) and T J k ) functions. In the case of the primitive model, using the electroneutrality condition C l p l q l= 0, (8) might be rewritten as In yi = - 2 ~ x p ~ S ~ ' ~ c *d, r~+( -,&kaik3 r2x) r ~ k
O
3
k
+
where Ui is a dimensionless average potential interaction energy, which could be attributed to an ion of species i
Similarly, the excess internal energy is given by the energy equation, which for temperature-dependent (solvent-averaged) pair potentials is given by
1
I
0.019 65
? J k ) = ( 4 x / k ) 1 - ~ * ~ , ( rsin ) r ( k r ) d r - Pqrql/tk2 ( 1 1)
- -
0.102 0.102 0.273 0.273 0.275 0.431 0.435 0.439 0.548 0.555 0.563 0.651 0.663 0.677
2 3
The Fourier transform of cl,(r) is then given by the formula
Similarly, it is found that yij(r) -cij(r) @ui,(r)as r and we thus define another short-ranged function by
0.101
I 2 3
0.970 0.970 0.971 0.945 0.947 0.950 0.979 0.992 1.008 1.093 1.134 1.190 1.353 1.480 1.672
-0.0967 -0.09 59 -0.09 52 -0.229 -0.218 -0.207 -0.259 -0.196 -0.136 -0.124 0.076 0.276 0.273 0.847 I .48
-0.0967 -0.0969 -0.0969 -0.229 -0.232 -0.232 -0.259 -0.277 -0.277 -0.124 -0.182 -0.187 0.273 0.1 12 0.092
where p = C,p,is the total particle density and N the total number of particles in the system. In the case of the PM, and setting &/aT = 0, ( 1 7 ) might be rewritten as E e x / N k T = (1 / P ) C P , U ,
(17*)
I
with U, given by (16). It should be pointed out that t in real electrolyte solutions is quite temperature dependent. It is easy, however, to correct Eex values obtained by (17*) if the temperature dependence of c is known.2' All the EeXresults given in this work apply for hypothetical electrolyte solutions with &/aT = 0. From F X / N k Tand the contact values of the radial distribution functions one might also calculate the osmotic coefficients (a) by the pressure equation, which for the primitive model is given byZ2 = 1
2*
1
+ - ~ ~ p , p , ~ ~ , ~ g ,+, (T ~E e, x, )/ N k T 3P
I
J
(18)
with E e x / N k Tgiven by (17*).
Results The HNC approximation was solved for a number of 1 :1 and 2: 1 electrolytes, with parameters corresponding to aqueous solutions at 25 O C . The 1:l electrolytes were specified by a dimensionless Bjerrum parameter B = P1q+q-1/4ma+- = 1.681 and the 2:l electrolytes by B = 3.400. The salt density might be specified by a dimensionless parameter ps* E pa+?/(u+ + u-), where u+ and v- are the stoichiometric coefficient of the cation and the anion, respectively. In the present calculations p * $ was varied between p * s = 0.000421 2 and 0.090 95 in the case of the 1 : l electrolytes and between p * s = 0.000421 2 and 0.124 in the case of the 2:l electrolytes. The ratio d+/d- = a++/a-- between the cation and the anion diameters was varied between 1 and 3 in the case of the 1 :1 electrolytes and between 0.5 and 2 in the case of the 2: 1 electrolytes. The HNC results are supplemented by a few MC results, some of which have previously been published in parts and some of which are new. Our MC program, which is based on the Metropolis procedurez3for the canonical ensemble, has been described e l ~ e w h e r e . ~ ~ , ~ ~ Thermodynamic Properties. H N C results for the activity coefficients, excess internal energies, and osmotic coefficients are given in Tables I and 11. The 1:l electrolyte H N C results for In y+ and In y- are plotted in Figure 1 as functions of P * , ' / ~ . The (21) Rasaiah, J . C.; Friedman, H. L. J . Chem. Pbys. 1969, 50, 3965. (22) Equation 3 in ref 4 was wrongly stated and should be replaced by
o= I
(18) See, e.g.: Stell, G. In Sfatistical Mechanics. Parr A : Equilibrium Techniques; Berne, B. J., Ed.; Plenum Press: New York, 1977. (19) Filon, L. N. G. Proc. R. SOC.Edinburgh 1928, A49, 38. (20) Mandel, F.; Bearman, R. J . J . Compur. Phys. 1971, 7, 637.
1
+ -E'"/NkT + ' ~ * [ g + + ( d + ) ( d + / a+)g~- - ( d - ) ( d - / a ) 3 + 2g+-(a)] 3
6
( 2 3 ) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 1087. (24) Sloth, P.; Ssrensen, T. S.; Jensen, J. B. J . Cbem. Soc., Faraday Trans. 2 1987, 83, 881.
Single-Ion Activity Coefficients
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2119
t"'
I"Y+
I
I
0.2
0.1
0.3
1.o
1.c
1.5
I-
1
!i- 1 ''W
1.0-
tI
\
\DHLL
Figure 1. Logarithm of single-ion activity coefficients for 1:l electrolyte solutions ( B = 1.68 1). Full lines, H N C results; filled circles, M C results for d,/d- = 1; squares, M C results for d+/d- = 3.
M C results from ref 4a are also indicated in the figures for comparison, as well as the limiting law of Debye and HuckelZS (DHLL). I t is seen that In y+ is quite dependent on the ratio d+/d-, at moderate to high concentrations, whereas In y - shows only a weak dependence on d+/d-. This indicates that the activity coefficient for the smallest ion in the system is rather unaffected by the specific interactions between co-ions (in the PM specified by a++ and a_-). The activity coefficients of the larger ion (here the cation), on the other hand, increase significantly as the ratio d+/d- = a++/a-- gets larger. ( 2 5 ) Debye, P.; Hiickel, E. Phys. Z.1923,24, 185.
Figure 2. Logarithm of single-ion activity coefficients for 2: 1 electrolyte solutions ( B = 3.400). H N C results for d+/d- = 0.5, 1, and 2 are indicated by full lines, dashed lines and dotdashed lines, respectively.
The HNC results for the single-ion activity coefficients of the 2:l electrolytes are plotted in Figure 2. It is seen, that In y+ does not depend strongly on d+/d- in contrast to In y-, which increases
2120
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990
TABLE 11: Values of the Excess Internal Energies, Activity Coefficients, and Osmotic Coefficients Obtained by the HNC Approximation for Aqueous 2 1 Electrolytes ( B = 3.400) P*< dJd-ECXINkT 0 In y+ In y.. 0.000421 2 0.5 0.345 0.905 -0.613 -0.175 0.344 0.904 -0.613 -0.177 1 -0.177 2 0.344 0.904 -0.613 0.805 -0.336 0.004 797 5 0.5 0.845 -1.427 1 0.802 0.834 -1.430 -0.366 -0.375 2 0.802 0.830 -1.431 -0.178 1.171 0.943 -1.983 0.019 65 0.5 I -0.350 1.160 0.878 -1.992 -0.401 1.154 2 0.86 1 -1.983 1.295 -2.152 0.536 1.438 0.046 24 0.5 1.060 -2.152 -0.062 1 1.409 1.010 -2.086 -0.229 1.393 2 2.464 -1.779 2.857 1.696 0.090 95 0.5 0.726 1.628 I 1.538 -1.774 0.225 1.600 2 1.396 -1.491 6.15 1.84 4.21 -1.04 0.124 0.5 1.57 2.08 -1.12 I 1.74 1.82 0.666 1.70 -0.547
significantly as the ratio d + / d - decreases. The In y+ curves for d + / d - = 0.5 and d + / d - = 1 are in fact almost coincident at the present scale. We have also performed a (rather large) number of MC runs for the most dense equal radii 2:l electrolyte studied in this work ( P * ~= 0.124). The total number of ions (N) in the MC cell was varied between N = 9 and N = 600, and at least a total number of 2 million configurations were generated in each run. It was found, however, that the variation of the thermodynamic properties with the size of the MC cell is rather complex in this case, making an exact extrapolation to the thermodynamic limit ( N m) somewhat difficult. This can readily be seen from Figure 3, where our MC results for E C X / N k Tare given as a function of 1 / N . In this figure, the corresponding grand canonical ensemble MC results of Valleau and Cohen26 are also indicated for comparison. It should be mentioned that the values of ECX/NkTobtained in each run did not vary significantly during the last 1 million configurations, indicating that the behavior shown in Figure 3 cannot be explained simply by statistical scatter. The M C results for the largest systems are, by the way, in very good agreement with the corresponding H N C result, but clearly more work has to be done in order to establish correctly a more accurate MC result in the limit N a. Radial Distribution Functions. Radial distribution functions for the 1:1 electrolytes with d+/d- = 1 have been published quite some years ago,*' and here we discuss primarily some of the most interesting features found for the nonequal radii 1: I electrolytes. In Figure 4 are shown the g+-(r) functions for p*s = 0.09095 and d + / d - = 1 and 3, obtained by the H N C approximation and by MC calculations. First it is noted that the agreement between the HNC and MC functions is excellent in these cases. It is seen that the g + - ( r ) function does not depend very much on the ratio d+/d-, except for the distinct minimum at r 2.3a+-found for d+/d- = 3. It can be mentioned that a similar minimum is found for the corresponding function in a similar uncharged hard-sphere system. At the lower densities, the g+-(r) functions are even more In Figure unaffected by the ratio d+/d-, for a given value of 5, g+ + ( r )functions calculated by the H N C approximation for B = 1.68 1 and d + / d - = 3 are shown for p * s = 0.019 65,0.046 24, and 0.090 95. For comparison we also show the corresponding MC result for the most dense system. The g+ + ( r ) function for p * s = 0.01965 is monotonic increasing-at least in the interval shown here. This implies that the repulsive electrostatic force dominates g , + ( r ) at short range at this concentration. As the density increases, oscillations begin to appear, and at the highest density ( P * ~= 0.09095) g + + ( r ) is changing into a decreasing function near contact. This indicates that the effective attractive
-
Sloth and Smensen
'1
1.90
1.85--
1.80-
*
*
"i..'
,
,
,
0,blO
,
I
,
O,b20
0,030 VN
Figure 3. MC results for Eea/NkTas a function of I / N for the 2:l electrolyte solution specified by B = 3.400, p * s = 0.124, and d+/d- = 1. 0 , canonical ensemble results; V,grand canonical ensemble results from ref 26. The value in parentheses corresponds to the slightly lower density; p*s 3 0.123.
-
(26) Valleau, J. P.; Cohen, L. K. J . Chem. Phys. 1980, 72, 5935. (27) Card, D. N.; Valleau, J. P. J . Chem. Phys. 1970, 52, 6232.
I
I
I
I
l-
1
2
r/a+--i
Figure 4. g + _ ( r ) functions for 1:l electrolytes with B = 1.681 and p * % = 0.090 95. Dashed line and squares, HNC and MC results for d+/d= I ; full line and filled circles, HNC and MC results for d+/d- = 3.
hard sphere (averaged) forces contribute significantly to the g++(r) functions at high densities.
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2121
Single-Ion Activity Coefficients
9, - ( r/a* - )
1.3
\\ b
1.2
\
\
1.1
1.o
-
r/a+ -4, I
1
2
Figure 5. g+t(r) functions for 1:l electrolytes with B = 1.681 and d+/d_ = 3. Short dashed line, HNC result for p*a = 0.01965; full line, HNC result for p l r = 0.04624; long dashed line, HNC result for p*s = 0.09095; 0 , MC results for = 0.09095.
Figure 6. g t - ( r ) functions for 2:l electrolytes with B = 3.400 and P * ~ = 0.124. Dashed line, full line, and filled circles indicates HNC results for d t / d - = 0.5, 1.0, and 2.0, respectively.
peculiar oscillating g--(r) function. All of the g--(r) functions for the 1:l electrolytes investigated here are increasing with distance near contact, showing only minor oscillations at high density. The g+-(r) functions for the most dense 2:l electrolytes ( P * ~ = 0.124) are shown in Figure 6. It is seen that these functions are quite independent of d+/d-, the most significant difference being the more pronounced oscillations which occur for d+/d- = 0.5. As the density decreases, the g+-(r) functions become more and more independent of d+/d-, as is the case for the 1:l electrolytes. Radial distribution functions between co-ions are shown in Figure 7 and 8 for p*s = 0.004 797 5 and 0.124. At the lower concentration, we see that these functions-as well as g+-(r)-are almost independent of d+/d-, except for the different “cutoff“ caused by the hard cores. At the high concentration, however, the g++(r)and, especially, the g J r ) functions show a significant dependence of d+/d-. All of the radial distribution functions shown in Figure 8, except g--(r) for d+/d- = 0.5, have a local maximum at r = 2.la+-. These maxima correspond roughly to configurations where two co-ions are placed on a line, separated by a counterion. The shape of the g--(r) function for p*s = 0.124 and d+/d- = 0.5 differs somewhat from this pattern, but is in fact very similar to the corresponding function in a neutral hard-sphere system, which is indicated by triangles in Figure 8 (bottom). In Figure 9, the H N C results for the g--(r) functions are shown for the three most dense, equal-radii, 2:1 electrolytes studied in this work. For comparison are also shown M C results for p * s = 0.09095 and 0.124. The H N C approximation is found to reproduce the MC function reasonably well for p*s = 0.124, whereas some deviations are found for p * s = 0.09095 at short distance. At the lowest density shown in this figure ( P * ~= 0.046 24), the g--(r) function is monotonic increasing at short range ( r / a + - < 2.2), implying that the repulsive electrostatic forces dominate whenever two anions are close to each other. For p*s = 0.090 95, the repulsive electrostatic and the (effective) attractive hard sphere forces between two anions seem to compete, which results in a rather
Application to KCI Solutions One might ask whether or not the primitive model can be used as a reasonable model for real electrolyte solutions. Here we recall that the PM in fact involves two successive approximations. First the McMillan-Mayer (MM) level14of description is applied: The solvent effects are included implicitly by the use of solvent-averaged pair potentials, Le., the potentials of mean forces between two solute species at infinite dilution, and the solvent-averaged ion-ion interactions are assumed to be functions of pair potentials only. Second, one approximates these solvent-averaged pair potentials by (1). Equation l yields the correct solvent-averaged pair-potential functions for large r. For small values of r, however, the PM is a quite poor approximation which does not describe the oscillatory behavior expected to occur for real electrolyte solutions. Recently, Kusalik and Patey investigated some extended models-which treat the solvent on a molecular level (the so-called Born-Oppenheimer (BO) level)-for a number of electrolyte solutions within the (R)HNC approximation.5 In these models the ions were modeled by hard charged spheres with radii estimated from X-ray electron density measurements.28 The solvent was represented by “water-like” polarizable polar hard spheres. In the case of a model KCI solution it was found that the application of solvent-averaged pair potentials (which may consistently be estimated for the model considered) yielded results close to those found by including the solvent on a molecular level. The activity coefficients obtained by the two different methods were in fact very similar for salt concentrations as high as 4 mol/dm3. This indicates that, in some cases, the application of solvent-averaged pair potentials might yield quite reasonable results for real systems-at least if the density is not too high. The next question to be asked is then: Could reasonable results for the (28) Morris, D.F. C. Srruct. Bonding 1968, 4, 63
2122
The Journal of Physical Chemistry, Vol. 94, No. 5 , 1990
f
Sloth and Serrensen
g++wa+.)
0.51
i.
2
0 0 0 0 0
I I
Figure 7. HNC results for radial distribution functions between co-ions for 2:l electrolytes with B = 3.400 and p*s = 0.004 797 5 . Dashed lines, d+/d_= 0.5; dotted lines, d+/d_= 2.0.
thermodynamic properties of a real electrolyte solution be obtained by approximating the solvent-averaged pair-potential functions by ( I ) ? Although the radial distribution functions obtained by
I
2
3 r/a+-
Figure 8. HNC results for radial distribution functions between co-ions for 2:1 electrolytes with B = 3.400 and p*s = 0.124. Dashed lines, d+/d= 0.5; full lines, d+/d- = 1; dotted lines, d+/d- = 2. Triangles indicate the ‘‘g..-(r)” function obtained by the HNC approximation for a hardsphere fluid specified by B = 0, p * s = 0.124,and d+/d_= 0.5.
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2123
Single-Ion Activity Coefficients
I
-0.5
t
t
t +
Figure 10. In y , vs square root of salt concentration (in mol/dm3) for KCI solutions. Triangles are the MM level results for the extended model of Kusalik and Patey (D = 933.' The full line is the corresponding PM results; the dashed line is the PM results for D = 78.5 and filled circles are the experimental values.
."
\
I
1
2
r/a, -1
the ratio d+/d- is close to unity. It should be mentioned that the results taken from ref 5 were obtained by application of the Kirkwood-Buff formalism,30which relates the derivative of In y , with respect to the salt density to structural properties. This approach has the disadvantage that one has to solve the HNC equation over a range of concentrations, which is not the case by application of (8), where the result for a given concentration might be obtained within a single calculation. From Figure 10 it is seen that the extended model yields activity coefficients that are very different from the corresponding primitive model results ( D = 93.5). Furthermore, we see that the primitive model results for D = 78.5 are in surprisingly good agreement with the experimental data for aqueous KCI solutions at 25 OC. Thus it seems as if the PM, in fact, offers a better description of the activity coefficients in real KCI solutions than the extended model. It should be mentioned, however, that the H N C approximation may fail to describe mixtures of ions and dipoles correctly: and this fact might explain the large difference between the extended model results and experimental data. Clearly, more work has to be done in order to obtain a reliable description of specific solvent effects in real electrolyte solutions. Finally, it should be noted that we did not try to adjust the constant solvent activity data obtained by the primitive model to the corresponding constant pressure data. This correction is small, however, at densities that are not too high. Indeed, one might expect the primitive model to be rather inadequate at densities where this correction becomes important. Acknowledgment. Financial support from the Danish Natural Science Research Council and Teknologistyrelsen (The National Agency for Technology) is gratefully acknowledged. ~~
(29) Hamer, W. J.; Wu, Y . C . J . Phys. Chem. Ref. Data 1972, I, 1047.
~~~
(30) Kirkwood, J. G . ; Buff, F. P.J . Chem. Phys. 1951, 19, 774.
