J . Phys. Chem. 1987, 91, 2978-2981
2978
For lower KR,1/yi2 increases; that is, y+ decreases but remains within physically reasonable values. The values of ?* predicted by the K equation become too low for small values of K R ,where the presence of the counterions is significantly important. I
I
I
’Q
I
1
i
15
lb
YR Figure 3. The reciprocal mean squared activity coefficient l / ~ plotted , ~ as a function of the reduced charge parameter y R for values of the reduced salt concentration parameter K Rof 1, 10, 50, and 100. The continuous curve corresponds to the variational solution of the KY equation discussed in this paper. The dashed curve is the result of the K
equation. Though not plotted in this figure, for very large values of p2 all three curves converge to the same values, represented by the Debye-Hiickel or Rice-Whitehead limit. In Figure 3 we show the behavior of I/yk2as a function of the reduced surface charge y R for several values of KR. We can see that for values of K Ras large as 100 the mean activity coefficient is essentially equal to one for all values of y R in the range shown.
Conclusions Starting form the expression of the electrostatic free energy associated with the Poisson-Boltzmann equation and using the surface charge density as the independent fixed variable, we obtained a new expression for the mean electrostatic coefficient of an electrolyte solution inside a cylindrical charged micropore (eq 9). In writing the Poisson-Boltzmann equation for this system, we have given particular attention to the presence of counterions inside the micropore and obtained an equation that contains two terms, the first one associated with the added salt ( K term) and the second one, the y term, associated with the counterion concentration and hence to the surface charge density. We have solved this so-called KY equation variationally, using the method of Arthurs and Robinson by minimizing a functional that corresponds to the electrostatic free energy. A comparison to a numerical solution of the differential equation showed the good accuracy of this approximate solution. The values obtained for the potential at the contact, c $ ~ and , for the mean ionic coefficient for the ion-surface interaction showed the importance of the screening introduced by the counterions. In a future report we shall present a numerical solution of the KY equation and discuss the effect on other important phenomena as the so-called ionic exclusion and the electrokinetic flow. Acknowledgment. We would like to acknowledge many helpful discussions with Pedro Colmen>res and the financial support of the Consejo de Desarrollo Cientifico (CDCH), of the Universidad de Los Andes (Grant No. C-235).
Solvation of Multiply Charged Ions: Predictions Using the Reference Interaction Site Method with Hypernetted Chain Closure Robert A. Kuharski and David Chandler* Department of Chemistry, University of California, Berkeley, Berkeley, California 94720 (Received: December 23, 1986)
The extended reference interaction site method using the hypernetted chain closure (HNC-REM) has been applied to study solvation structure and thermodynamics of di- and trivalent ions in aqueous solution. The theory predicts incorrect values for numbers of solvent neighbors surrounding highly charged ions, and this inaccuracy leads to incorrect predictions for the activation energy for adiabatic electron transfer. We consider Fez+ Fe3+,and in this case, the prediction is about 20 kcal/mol too high. Errors of this type as well as violations of stoichiometric conditions reveal the need for an improved closure in the RISM theory for ionic solvation.
I. Introduction The reference interaction site method (RISM) integral equation in its extended form (i.e., with the H N C closure)z has been applied to the study of liquid water structure3 and ionic solution^^*^ The results obtained in this way by Rossky and co-workers represent an important development since these calculations demonstrate the possibility of understanding aqueous solutions and ionic solvation with a computationally convenient theory. Motivated and (1) Chandler, D. In The Liquid State of Matter: Fluids, simple Complex; North-Holland: Amsterdam, 1982. (2) Hirata, F.; Rossky, P. J. Chem. Phys. Lett. 1981, 83, 329. (3) Pettitt, B. M.; Rossky, P. J. J . Chem. Phys. 1982, 77, 1451. (4) (a) Pettitt, B. M,; Rossky, p. J , J , Chem, phys. 1986, 84, 5836, (b) Hirata, F.; Rossky, P. J.; Pettitt, B. M. J. Chem. phys. 1983, 78, 4133. ( 5 ) Chiles, R. A.; Rossky, P. J. J . Am. Chem. SOC.1984, 106, 6867.
0022-3654/87/2091-2978$01.50/0
by this development, we have examined the predictions of the HNC-RISM heory for the aqueous solvation 2 multiply charged ions, and this paper reports on our findings. Most significant, we find that the HNC-RISM equation overestimates the changes in solvation structure which occur due to changes in ionic charge. We examine the ramifications of this overestimate and suggest the necessity of improving the theory when applied to such systems. At the outset, it should be noted that problems with the HNC-RISM theory are already apparent in the literature. For instance, Pettitt and Rossby’s4 calculations for numbers of neighbors of water molecules surrounding monovalent ions reveal an inconsistency in that the numbers calculated from gaH(r)are different than those determined with gao(r). Here, gas(r)refers to the HNC-RISM radial distribution function between ion cy and atom s of a solvent water molecule. In the solvation of Na+ ions, Pettitt and Rossky compute 4.28 neighbors from the distribution 0 1987 American Chemical Society
1rhe Journal of Physical Chemistry, Vol. 91, No. 11, 1987 2919
Solvation of Multiply Charged Ions of oxygen atoms surrounding a Na+ ion, but they find 6.97 neighboring water molecules when considering the distribution of protons. This disparity can be viewed as representing a failure of the theory to properly adhere in the small length scale regime to the constraints of stoichiometry. (In the long wavelength limit, all RISM theories by construction obey the constraints of stoichiometry.) Another problem is hinted at in the Chiles and Rossky calcul a t i o n ~of~the solvation contribution to the free energy of activation for the SN2 reaction CI- + CH3CI C1CH3 + C1-. In this case, computer simulations6have determined that the free energy barrier in aqueous solution is about 26 kcal/mol while for very nearly the same model system the HNC-RISM theory predicts 31 kcal/mol. This type of overestimate obtained from the theory is, we will see, significantly magnified if one considers multiply charged ionic species, and the large overestimates are directly related to incorrect predictions of solvation numbers. In our own work (see section 11), we have found a curious dependence of the pure solvent water HNC-RISM pair distribution functions on variations of the site-site potentials in physically inaccessible regions. The dependence is far greater than might be anticipated from earlier studies of the effects of auxiliary sites in the RISM theory.' To examine the predictions of this theory for ionic solvation, we consider a model constructed to mimic Fez+ and Fe3+ ions in water. The necessary equations and site-site potential are discussed in section 11. The results of the theory are presented in section 111.
-
11. RISM Equations and Models
We will consider the solvation of one ion in water or two ions separated by a distance R in water. In the latter case, of relevance to the study of electron transfer, the solute species can be viewed as a rigid charged diatomic species at infinite dilution in water. For this situation, the Ornstein-Zernike-like equation introduced by Chandler and Andersen* is
ph,,(k) = D , , ( k ) Z,,4k) 2 d k )
(1)
Y.S'
where a = 1 or 2 refers to the two different ions sin ( k R ) q 2 ( k )= kR B , , ( k ) = 0 2 , ( k ) = 1, p is the bulk number density of the solvent, the Roman subscripts refer to the interaction sites of the solvent, and ~ , , , ( kis) the Fourier transform of the site density pair correlation function of the pure solvent. The circumflexes refer to the Fourier space representation. The Chandler-Andersen equation contains two sets of unknown functions: the pair distribution function ges(r) = 1 hus(r),and the direct correlation function, c,,(r). To close the equation, an additional relationship is required. The H N C closure is C,&) = --@,Ad + ha&) - In [1 + h u m 1 (3) where u,,(r) is the pair potential between ion a and solvent site s, and as always, p' = kBT. After specifying the potential and the pure solvent correlation functions, eq 1 and 3 can be solved self-consistentl y. Once the cmS(r)functions are determined in this way, the solvation contribution to the solute chemical potentials, pLsolv, can be computed from
+
0 ' 0.
I
2.
4.
6.
8.
10.
Figure 1. Radial distribution functions for Fez+and Fe3+in liquid water. The solid and dashed lines are the ion-oxygen and ion-hydrogen distributions, respectively. The numbers under the arrows indicate the height of the main peak in the ion-oxygen distributions.
The pertinent equations for studying the solvation of just one ion are eq 1, 3, and 4, where the indices a and y refer to only one site. The pure solvent correlation functions gSS,(k)are needed in the implementation of these equations. Simulations or experiments could provide this information. We, however, have obtained these functions by solving the corresponding HNC-RISM equations for liquid water at 25 OC and 1 g/cm3 with the intermolecular potentials of the SPC form. In particular, we have augmented the SPC model of Berendsen et al.1° in precisely the same way as done by Pettitt and R ~ s s k y .The ~ alteration adds c(a/r)lZto the 0-H intermolecular site-site potential of the SPC model, where t = 0.8 kcal/mol and a = 1.6A. An addition of this type is required to ensure stability of the model since without it the total waterwater pair potentials possess configurations with --m potential energy. Since these regions are surrounded by large potential maxima, it would seem that any reasonable choice of potential added to create a stable model would yield the same correlation functions as any other reasonable choice. However, we have found that the position of the first peak in the HNC-RISM gHO(r)is unduly sensitive to that choice. For instance, changing a from 1.6 to 1 8, causes virtually no change in the intermolecular potential, but it alters the location of the first peak in the HNCRISM go&) by 0.5 8,. In face of this difficulty, we have adopted the potential suggested by Pettitt and Rossky and described above. For the interactions between ions and water, we have used
u,,(r) = 6 , , ~ ( ~ / +r z,zse2/r )~
(5)
where E = 12.5 kcal/mol, L = 2 A, 60, is 1 when s refers to oxygen and zero otherwise, z,e is the charge of the ath ion, and z,e is the partial charge of site s in the SPC model of water. The parameters in this model predict a water-ion minimum energy distance of 1.96 and 1.85 A for z , = 2 and 3, respectively. These distances can be compared with the corresponding lengths of 1.97 and 1.85 A for Fez+ and Fe3+,respectively, suggested by Newton's calculations on ferric and ferrous water clusters." Note that the model was constructed to give the lengths as described here, but little attention was paid to the absolute value of this energy. In a more realistic model, one must consider effects such as the polarization of water molecules and the concomitant many-body interactions. These effects undoubtedly play a significant role in the absolute value of the energy.
This formula results from the standard coupling parameter integration when applied to the HNC-RISM t h e ~ r y . ~
111. Results Equations 1 and 2 were solved self-consistently with the potential models described in the previous section. The numerical method
( 6 ) Chandrasekhar, J.; Smith, S. F.; Jorgensen, W. L. J . Am. Chem. SOC. 1984, 106, 3049. (7) H s u , C. S.; Chandler, D.; Lowden, L. J. Chem. Phys. 1976, 14, 213. Cummings, P. T.; Gray, C. G.; Sullivan, D. J . Phys. A 1981, 14, 1483. (8) Chandler, D.; Andersen, H. C. J . Chem. Phys. 1972, 57, 1930.
(9) Singer, S. J.; Chandler, D. Mol. Phys. 1985, 55, 621 (10) Berendsen, H. J. C.; Postma, J. P. M.; Von Gunsteren, W. F.; Hermans, J. Infermolecular Forces; Reidel: Hingham, MA, 1981. ( 1 1) Jafri, J. A.; Logan, J.; Newton, M. D. Isr. J . Chem. 1980, 19, 340.
The Journal of Physical Chemistry, Vol. 91, No. 11, 1987
2980
TABLE I: Solvation of a Single Ion in WateP coord no. ionic charge z ,
Ob
H'
-PPSOl"
- ~ A E ~
2 3 4
5.2 1.4 9.7 12.1 14.9
5.5 6.4 1.5 9.0 11.1
607.4 1460.8 27 14.9 4390.6 6505.1
1220 2901 5391 8129 12950
5 6
'All numbers computed from HNC-RISM integral equations with potential described in text. bComputed from ion-oxygen radial distribution.
