J . Phys. Chem. 1990, 94, 21 13-21 16 is also achieved under conditions of relatively lower salt concentration at small surface charge density. In this case, however, it is more likely to be due to a significant correlation effect due to co-ions.
2113
At sufficiently high surface charge density and salt concentration, and in the presence of divalent ions, there is a charge reversal in the double layer with a possibly attractive net osmotic pressure even for a single double layer.
Dynamlcal Effects of Ions on a Polar Solvent S.-B. Zhu, J. Lee, and G . W. Robinson* Picosecond and Quantum Radiation Laboratory, Texas Tech University, P.O. Box 2460, Lubbock, Texas 79409 (Received: March 31, 1989; In Final Form: September 25, 1989)
The dynamical behavior of ion pairs in a polar diatomic solvent is investigated by using molecular dynamics (MD) computer simulations. Two methods are employed. One is a full-scale MD calculation without approximations, while the other uses the clamping approximation (CA). A comparison of the results from these two methods gives information about the dynamical effects. For example, it is found that for LiF the cation attracts more solvent molecules than the anion does. This effect is enhanced by the mass difference of Li and F on the dynamics. This result is in line with experimental observations. The relaxation features of the rotational and translational motions of the neighboring solvent molecules are quite different. The former exhibits a high-frequency oscillation induced by the high-frequency intramolecular vibrations of the ion pair, while the latter decays as a damped oscillator with a much lower frequency. A nonexponential response of the molecular reorientation of the polar solvent is observed.
Introduction Extensive molecular dynamics (MD) simulations’-5 have been focused on variations of the local structure and properties of the solute under the influence of the solvent. However, a knowledge of the effects of a solute molecule on the neighboring solvent is also required for a complete understanding of intricate molecular mechanisms in many condensed-phase phenomena, including chemical reactions. Therefore, in this paper, we continue this type of investigation but shift our attention to a less discussed, nevertheless equally important, aspect: dynamical effects of ions on a polar solvent. Classical Brownian motion theory deals with massive solute particles. In this limiting case, one can simplify the problem of reaction dynamics by using a “clamping approximation” (CA).6 However, this simple method cannot be adequately extended to certain problems that depend on ultrafast dynamical properties of the system.’ In certain circumstances, the missing dynamical effects of the solute may lead to totally different conclusions than are obtained by using a more exact treatment. To confirm this point further, in this paper we will compare the results obtained from a full-scale molecular dynamics (MD) calculation with those obtained from CA. Method The ensemble under study includes 107 diatomic molecules representing a typical polar solvent plus one LiF molecule capable of dissociating into a Li+ ion and a F ion. Conventional periodic boundary conditions with an interaction cutoff distance of 8.1 13 A are used to extend the system. The diatomic molecule is modeled as a two-point entity, each point being taken to be a Lennard-Jones particle. For definiteness, one of these particles is taken to be nitrogen, the other oxygen. The nitrogen atom is assigned a small positive charge, while the oxygen atom is neg( I ) Geiger, A.; Rahman, A.; Stillinger, F. H. J . Chem. Phys. 1979, 70, 263. (2) Rao, M.; Berne, B. J. J . Phys. Chem. 1980, 85, 1498. (3) Impey, R. W.; Madden, P. A.; McDonald, I. R. J . Phys. Chem. 1983, 87, 5071. (4) Goldman, S.;Back, P. J . Chem. Phys. 1986, 84, 2761. ( 5 ) Zhu, S.-B.;Robinson, G. W. J . Chem. Phys. 1989, 90, 7127. (6) See, e.g., ref 4. (7) Zhu, S.-B.; Lee, J.; Robinson, G.W. Phys. Reu. 1989, ,439, 5985.
0022-3654/90/2094-2113$02.50/0
TABLE I: Molecular Parameters’ ULir A 2.948 qLi, CGSE
A UN, A OF,
A
00, CUlkB,
K eplks, K tN/kB, K fo/kB, K
3.123 3.394 3.131 39.3 112 91.5 113
qF, CGSE 9 ~ CGSE , 90, CGSE k,,
dyn/cm
re, A
4.8 X lo-’’ -4.8 X lo-’’ 1.39 X IO-’’ -1.39 X IO-” 1.594 X lo6 1.15
’The Lennard-Jonesparameters for different particles are estimated
by using the usual Lorentz-Berthelot combining rules.
atively charged to reproduce a dipole moment of 0.16 D. We write the solvent-solvent intermolecular potential as
where rlikjdenote the distance between atom 1 of the ith solvent molecule and atom k of the j t h solvent molecule; qIi and qk, represent the point charges carried by atoms 1 and k of the diatomic solvent. The solvent molecules are not rigid, the two atoms being linked through a harmonic potential
Both solute atoms are coupled to the atoms of the solvent molecules through Lennard-Jones and Coulomb interactions
with 1 representing Li or F and k representing the solvent atoms. The lithium atom and the fluorine atom of the solute are attracted to each other through the long-range electrostatic interaction and repel each other according to a short-range power law potential, ULi, = Ax-“
+ qLi+qF-x-’
(4)
This potential is often referred to as the “Pauling form”.* There (8) Pauling, L. J . Am. Chem. SOC.1927, 49, 765.
0 1990 American Chemical Society
2114
Zhu et al.
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990
L
!I-
i
A
N
m8
i
y i !
1
I
I
3700 4800 Distance
I
I
I
5900
I
(X 0.001A)
I
7000
I
Figure 1. Li-0 radial distribution functions (RDFs). Curve represents results from the full-scale MD; symbols are from the clamping approximation. Same for Figures 2-4. I
b-
Figure 3. Li-N radial distribution functions.
:;!lLJ C'
iL L 11
I
':
LB -
B
0,
.
..
'e
...e**
0 .
LI I _ L
i
I
f
1
3700 4800 Distance
1
I
5900
I
(X U . 0 U l A )
I
7000
I
3700 4800 Distance
Figure 2. F-0 radial distribution functions.
are many ways of determining the parameters in eq 4 from a variety of experimental data.9 Here we simply fix n = 10 and adjust A to match the minimum in the potential energy at the LiF gas phase equilibrium bond length x,, that is A = (qLi+)*x,9/10
(5)
where x, = 1.55 A.Io This gives too large a dissociation energy by a factor of about 1.5. The discrepancy does not matter for our current purpose since we are not intending to reproduce any real experiment. All the parameters used in this simulation are contained in Table I . The MD is performed at a temperature of 120 K with the length of the cubic sample being 16.23 A. This is equivalent to a liquid density of 1.269 g/cm3. The standard leap-frog algorithm" is employed to integrate the equations of motion with a time step equal to 1.3 fs. A total of 1.2 X I O 5 time steps are performed for collecting data. For comparison, we also carry out an approximate molecular dynamics calculation on the same ensemble using CA. Four sets of such simulations are performed, fixing the Li-F separations at 1.68, 1.943, 2.206, and 2.469 8.Each simulation lasts 58.5 ps (4.5 X IO4 time steps) to obtain statistically reasonable results. To compare the results with those obtained from the full-scale M D calculation, we average the four sets of data according to the probability of each configuration. Static Properties Plotted in Figures 1-4 are solute-solvent atom-atom radial distribution functions. The symbols represent the results from the CA calculations averaged over the four chosen configurations with weights 0.2322, 0.2368, 0.4815, and 0.0495 evaluated from the MD simulation. Large deviations may sometimes occur at
P.Ado. Solid State Phys. 1964, 16, 1. ( I O ) Brewer, L.; Brackett, E. Chem. Reo. 1961, 66, 425. ( 1 1 ) Verlet, L. Phys. Reu. 1967, 159, 98.
(9) Tosi, M .
5900
(X 0.001A) Figure 4. F-N radial distribution functions.
I
7000
TABLE 11: Shell Sizes and Averarre Numbers of Particles Contained' shell boundaries no. of 0 atoms shell no. rmln,8, rmax, 8, MD CA 1
0 0
2
5.01 3.88
3
shell no.
-
5.01 3.88 6.76 6.64
-
12.06 5.30 27.41 21.34
11.63 5.72 27.91 19.57
6.64
8.14
29.44
34.08
shell boundaries rmtm. 8,
1
0 0
2
4.76 4.26 6.14
3
r,... 8, 4.76 4.26 6.14 7.01 8.14
-
-
no. of N atoms MD CA 11.84 7.83 5.51 28.81 40.09
11.64 7.44 7.00 30.08 38.60
a First entries are for the cationic shells, and second entries are for the anionic shells.
short distances, and many of the subsidiary peaks and valleys are poorly reproduced by the CA. An immediately observable feature in these figures is that there exist clearly discernible solvent shells around the positive and negative ions. Such shells, especially at large distances, become a bit ambiguous for polar molecules in nonpolar liquids5 It is seen from Table I1 that the number of nearest-neighbor solvent molecules surrounding the positive ion considerably exceed those surrounding the negative ion. The smaller Li' ion (refer to Table I) attracts more solvent particles and forms a larger cluster shell (see Table 11). This must be attributed in part to a different partitioning of the Lennard-Jones and electrostatic interactions for the cation and anion species. However, this difference, though still prominent, becomes a little less if the CA method is employed. Evidently, the dynamical motion of the ions can be one of the causes leading to this deviation. The lighter lithium ion moves so rapidly that the surrounding solvent molecules
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2115
Dynamical Effects of Ions on a Polar Solvent
I
I
1040
I
I
2080
I
I
3120
I
I
I
4160
Time (X 0. lfs) Figure 5. Center-of-mass VACFs of the neighboring solvent molecules. Solid curve: longitudinal component; dotted (darker) curve: transverse components.
cannot fully keep up with it. As a result, the cationic shell broadens. We did not observe such phenomena in the study of a polar molecule in a nonpolar liquid, even though the dynamic effects of the solute were considered. These results can be used to explain the electrochemical phenomenon where the effective radii of low-mass hydrated ions in aqueous solution increase dramatically from their ionic radii in crystals.I2 As a direct consequence of this, membranes of living cells are generally more permeable to heavy ions than to light ions. In fact, many important physiological mechanisms are based on this phenomenon.12 Asymmetric concentration distributions were observed long agoI3 in experiments on hydrated ions and have been confirmed more recently1“-” by the use of modern time-resolved techniques. Smaller ions bind more water than large ions, and cations more water than anions. In comparison with these experimental observations, our simulation results are quite encouraging and are able to provide a microscopic picture of this phenomenon.
Dynamical Properties In this section, we discuss the dynamical behavior of “neighboring solvent” molecules, defined as those whose centerof-mass distances to either the cation or anion are less than 3.824 A. According to this definition, the average number of neighboring polar molecules is 9.85 from the M D calculation and 9.68 from CA. Figure 5 shows that the solvent center-of-mass velocity autocorrelation functions (VACFs) for the component along the reaction coordinate (Li+-F direction) and for those transverse to the reaction coordinate are almost identical. Both types of VACF decay to zero within 100 fs. They demonstrate pronounced negative portions, indicating strong caging effects from the neighboring solvent. The similarity between the components parallel and perpendicular to the reaction coordinate indicates a homogeneous environment for the center-of-mass translations of the solvent molecules. Contrary to what is found for the translational VACF results, a highly oscillatory velocity autocorrelation function is evident for the angular velocity (AVACF). This is shown in Figure 6. The oscillations mainly follow the rapid intramolecular vibrations of the ion pair. Thus, in this case, the solvent molecules see an inhomogeneous environment. Such high-frequency oscillations in correlation functions were previously observed in MD calculations of molecular liquids in the presence of a strong, one-dimensional alternating external force field.’* The oscillations are missing for nonpolar solvent molecules surrounding a polar ~ o l u t e . ~ (12) See, e.g., Moore, W. J. Physical Chemistry, 4th ed.; Prentice-Hall: Englewood Cliffs, NJ, 1972; Chapter 10. (13) Bockris, J. O’M. Q. Reo. London 1949, 3, 173. (14) Struis, R. P. W. J.; de Bleijser, J.; Leyte, J. C. J . Phys. Chem. 1987, 91, 1639. (15) Struis, R. P. W. J.; de Bleijser, J.; Leyte, J. C. J. Phys. Chem. 1987, 91, 6309. (1 6 ) Lee, J. J. Am. Chem. SOC.1989, 1 I I, 427. (17) Lee, J. J . Phys. Chem., in press. (18) Zhu, S.-B.; Lee, J.; Robinson, G. W. Phys. Reu. 1988, ,438, 5810.
1040
2080
3120
4160
Time (X 0. lfs) Figure 6. AVACF of the neighboring solvent molecules.
I
I
1040
I
I
2080
I
I
3120
I
I
I
4160
Time iX 0. lfs) Figure 7. O A C F of the neighboring solvent molecules.
It seems that the ion pair plays totally different roles in affecting the translational and rotational motions of the neighboring solvent. Reorientational motions of molecular liquids can be measured in dielectric dispersion experimet~ts.l~*~ They are more important than velocity changes in the solvation response. Figure 7 displays this reorientational behavior. The relaxation is obviously nonexponential. This reflects a disagreement with the prediction of continuum models, which use the concept of a homogeneous dielectric medium, viewing the solute as a spherical cavity containing a point charge.21 It is not surprising that none of these dynamical effects can be reproduced by CA, which predicts nearly nondecaying correlation functions of the local solvent through neglect of the dynamical motion of the ions. It would certainly be informative to compare the exact M D results displayed in Figures 5-7 with results obtained by more advanced approximate methods such as MTGLEa2* This will be a subject for future work.
Conclusions By comparing the results obtained from MD simulations, with and without CA, we have investigated the dynamical behavior of a polar solvent in the presence of an ion pair. It is found that the cation shell contains more solvent particles than the anion shell. This effect is mainly caused by the interaction characteristics, such as ionic size, but it also arises partly because of the rapid dynamical motions of the smaller ion. This result is in agreement with experimental observations in aqueous salt solutions, which show that the hydration number is larger for a cation than for an ani~n.’~-~’ It was also found that ion pairs play different roles in affecting translational and rotational motions of inner-shell solvent molecules. For rotations, but not translations, the polar solvent (19) 257. (20) (21) (22)
Bagchi, B.; Oxtoby, D. W.; Fleming, G. R. Chem. Phys. 1984, 86, van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1985, 89, 4181. Maroncelli, M.; Fleming, G. R. J . Chem. Phys. 1988, 89, 5044. Adelman, S. A. Adu. Chem. Phys. 1983, 53, 61.
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
