J . Phys. Chem. 1990, 94, 6483-6488
dology behind the simulations are directed toward generating heat bath models for stochastic dynamics simulations of electronically excited species trapped in helium bubbles. The present pressure-dependent simulations provide many insights. The behavior of classical high-pressure He in the condensed phase is as may be expected. A key role is played by high-pressure effects in causing a liquidlike structure without the liquidlike binding (reflected in the internal energy of binding being an excess positive quantity). As pressure is increased, quantum effects become less (e.g., at 10-14-GPa cases reported here), and possibly, they become relatively nonnegligible for lower pressures (0.3-1.4 GPa). The low-pressure VACF's indicate gaslike behavior with a long hydrodynamic positive tail whereas the high-pressure ones show the oscillatory structure characterizing a liquidlike condensedphase correlation. The RDF's lead to similar conclusions. It is seen that the MTGLE heat bath parameters are quite sensitive to (and exponentially dependent on) the pressure. The sensitivity of the frequency moment parameters for He(%) and He(%), in comparison, reveals the change in the condensed-phase dynamics occurring as a result of a sudden electronic transition to which the solvation structure has not adjusted. Electronic excitation leads to the sampling of lower frequency components of the solvating medium as reflected in the slower decay of dy-
6483
namical quantities upon electronic excitation. It is also found that an increase of pressure may introduce less variations in the translational dynamics of the excited He(%) than it does in the ground state, He('S). The present heat bath models can be employed in stochastic dynamics calculations of electronically inelastic collision processes23within the high-pressure 4He matrix. The additional input needed will be the friction kernels, and a sophisticated semiempirical mode-coupling scheme, as employed by Gaskell and co-workers, is recommended for this purpose. However, a simple Gaussian model friction based on the present parameters has been employed in our preliminary stochastic dynamics studies.24 Acknowledgment. This work was supported by the Air Force Astronautics Laboratory at Edwards Air Force Base under Contract No. F04611-86-C-0068. It is a pleasure to thank Dr. Michael J. Redmon for his interest and encouragement. (23) (a) Garrett, B. C.; Swaminathan, P. K.; Murthy, C. S.; Redmon, M. J. J . Chem. Phys. 1987,87, 3207. (b) Swaminathan, P. K.; Garrett, B. C.; Murthy, C. S . J. Chem. Phys. 1988, 88, 2822. (24) Swaminathan, P. K.; Redmon, M. J.; Murthy, C. S.;Garrett, B. C.; Natanson, G. A. Dynamics of Electronically Excited Species in Gaseous and Condensed Phase: AL-TR-89-051, Dec 1989.
Solvent Polarization and Hydration of the Chlorine Anion Michiel Sprik, Department of Physics, Rijksuniversiteit te Utrecht, Postbus 80.000, 3508 TA, Utrecht, The Netherlands
Michael L. Klein,* and Kyoko Watanabet Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 (Received: October 23, 1989; In Final Form: February 20, 1990)
The effect of solvent polarization around a C1- ion in water is investigated by molecular dynamics simulation. Two different solvent models that each take account of the collective polarization are compared to the popular TIP4P potential. The first is an effective pair potential which incorporates polarization in a mean-field approximation. The second is a polarizable model which treats the many-body induction as a separate dynamical degree of freedom for each molecule. The hydration number found in the TIP4P solvent is reduced in both polarization models. Although the ion-solvent radial distribution functions for the polarizable models are in reasonable agreement with neutron scattering data, the time scale of the solvent relaxation processes in both models is too slow. This finding is in marked contrast with the case of TIP4P, where the time scales for such processes are too fast. The enhanced rigidity of the polarization models is a consequence of the large self-polarization correction.
I. Introduction One of the major applications of the molecular dynamics method is the study of strongly associative polar liquids among which water is the prime example. Many of the complex structural features of these systems such as hydrogen bonding and also dynamical behavior can be modeled by relatively simple molecular interaction potential^.'-^ The polar molecule is commonly represented as a rigid frame consisting of sites carrying fractional charges for the electrostatic interactions and sites where the short-range interactions of the atom-atom type are centered. Usually, the values of the fractional charges are chosen to give an effective dipole moment which is enhanced over the gas-phase value in order to account for the collective polarization effects in the bulk liquid. In this respect these models can, therefore, be considered as incorporating polarization using mean-field theory. Recently, it was realized4 that in a fully consistent treatment the self-interaction energy corresponding to the induced 'Present address: BP Research Centre, Sunbury-on-Thames, Middlesex TW16 7LN, United Kingdom.
0022-3654/90/2094-6483$02.50/0
polarization must also be included in the expression of the full internal energy of the liquid phase. As a consequence, the effective dipole moment is increased even more and now, for these so-called extended interaction site models, the equilibrium dielectric properties of water are also in fairly good agreement with experiment.5 The mean-field models are, therefore, quite successful in describing the properties of pure and homogeneous states. However, in situations where the distribution of molecules is inhomogeneous, e.g., at liquid-gas interfaces or near solutes, the mean-field approximation can be expected to be less accurate. The reason is (1) Stillinger, F.H.; Rahman, A. J. Chem. Phys. 1974,60, 1545. Skipper, N. T.; Neilson, G. W. J . Phys.: Condens. M a t e r 1989, 1, 4141. (2) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, H. J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; D r
331.
(3) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (4) Berendsen, H. J. C.; Grigera, J.; Straatsma, T. P. J . Phys. Chem. 1987, 91, 6269. ( 5 ) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157.
0 1990 American Chemical Society
6484
The Journal of Physical Chemistry, Vol. 94, No. 16, 199(3
that the average dipole moment in the boundary layer or the solvation shell is not necessarily identical with the value for the bulk. Indeed, already in the elementary case of the CI- anion solvated in water the agreement with experiment is less than satisfying. Most solvent models predict a coordination number of 7 or more,611while the experimental value is around 6.12 Also, the residence times of water molecules in the first solvation shell are found to be considerably longer than inferred from N M R measurements (see ref 8 and references therein). It should be pointed out that some potential models fitted, in part, to ab initio quantum calculations do yield hydration numbers around 6.13 A possible explanation for the overestimation of the coordination by the commonly used mean-field models6-11may be connected with the approximation of an average effective dipole moment for all solvent molecules that does not allow the induced polarization of the first solvation shell to adequately adapt to the local situation. One may speculate that the repulsive interaction between the near-parallel dipole moments, all pointing toward the anion, all lead to a depolarization and, therefore, to a reduction of the binding to the ion. Alternatively, the induced polarization may remain closer to the bulk value, but one or more molecules may be expelled from the solvation shell. In the latter case, the solvation may, perhaps, still be described by a mean-field model, but with a modified effective dipole moment. To further investigate these questions, in the present paper, we study the solvation of a single CI- ion in water using a recently developed polarizable model for waterI4 that includes the induced polarization as an explicit degree of freedom. We have also employed a related extended mean-field modelS in order to investigate to what extent an enhanced effective dipole moment and improved dielectric properties are able to account for the interactions in the solvation shell. The next section contains a review of the two models for pure water and some new results for the pressure and the dynamical behavior of the polarizable model.14 In section 3 the solvation of the CI-ion is treated, and we conclude in section 4 with a discussion. 11. Pure Water A . Polarizable Water. The development of interaction po-
tentials for water that model the induction of electronic polarization has a long history with varying degrees of s u c c e ~ s , I ~and -~~ recently there has been renewed activity in this field.I4Jg2’ Several different schemes have been proposed to handle polarization effects, but we will restrict ourselves to the model which we introduced in a recent p~blication.’~ The interested reader should consult the cited references for information about alternative approaches. The model that we employ is semiempirical and is based on (6) Impey, R. W.; Madden, P. A.; McDonald, I. R. J . Phys. Chem. 1983, 87, 5071. (7) Chandrasekhar. J.; Spellmeyer. D. C.; Jorgensen, W . L. J . Am. Chem. Soc. 1984, 106, 903. (8) Bopp, P. In The Physical Chemistry of Aqueous Ionic Solutions; Bellissent-Funel, M.-C., Neilson. G.W., Eds.; Reidel: Dordrecht, 1987; p 217. (9) Reddy, M . R.; Berkowitz, M . J . Chem. Phys. 1988, 88, 7104. (IO) Straatsma, T. P.; Berendsen, H. J. C. J . Chem. Phys. 1988,89,5876. ( I 1 ) Mezei, M.; Beveridge, D. L. J . Chem. Phys. 1981, 74, 6902. (12) Powell, D. H.; Barnes, A . C.; Enderby, J. E.; Neilson, G. W.; Salmon, P . S . Faraday Discuss. Chem. Soc. 1988,85, 137. Enderby, J. E.; Cummings, S.; Herdman, G.J.; Neilson, G . W.; Salmon, P. S. Skipper, N . T. J . Phys. Chem. 1987, 91, 5851. (13) Bounds, D. Mol. Phys. 1985, 54, 1335. Clementi, E.; Barsotti, R. Chem. Phys. Lett. 1978, 59, 21. (14) Sprik. M.; Klein, M . L. J . Chem. Phys. 1988, 89, 7556. (15) Stillinger, F. H.; David, C. W . J . Chem. Phys. 1978, 69, 1473. (16) Barnes, P.; Finney, J. L.; Nicholas, J. D.; Quinn, J . E. Nature 1979. 282, 459. (17) Levesque, D.; Weiss, J.-J.; Patey, G.N. Mol. Phys. 1984, 51, 333. Caillol, J . M.; Levesque, D.; Weiss, J.-J.; Kusalik, P. G.; Patey, G.N . Mol. Phys. 1985, 55, 65. (18) Rullmann, J. A. C.; van Duijnen, P. Th. Mol. Phys. 1988, 63, 451. (19) Ahlstrbm, P.; Wallqvist, A.; Engstrom, S.; Jonsson, B. Mol. Phys. 1989. 68, 563. Wallqvist, A,; Ahlstrom, P.; Karlstrom, G.J . Phys. Chem. 1990, 94, 1649. (20) Kuwajima, S.; Warshel, A . J . Phys. Chem. 1990, 94, 460. (2 I ) Modern Techniques in Computational Chemistry; MOTECC-89; Clementi. E , , Ed.: Escom: Leiden, 1989; p 414.
Sprik et al. a modification of the TIP4P potentiaL3 The geometry of the TIP4P monomer is retained, but the positive Q H and negative QM charges are scaled back to Q H = 0.443 e and QM = -0.886 e, respectively, in order to reproduce the experimental gas-phase dipole moment of 1.85 D (see Table I). In addition, the negative charge site is resolved into four separate sites (qi)arranged in a rigid tetrahedral configuration. Two sites are placed approximately along the OH bonds and the other two above and below the molecular plane, pointing in the “lone-pair” directions. The dimensions of the tetrahedron are kept small: the distance between two vertices is 0.283 A. The positions of the four charge sites are held fixed in the molecular frame. However, the values of the charges are allowed to fluctuate, subject to the constraint of a constant sum equal to the negative charge QM given above. The combination of the four charges can form an induced dipole of arbitrary magnitude and direction residing on the site of the negative charge. Mathematical details of the model can be found in ref 14. (An unfortunate, but obvious, error of signs in the orthogonal transformation matrix in eq 13 of ref 14 should be noted.) The fluctuating charges on each molecule are treated as regular classical dynamical degrees of freedom similar to the position coordinates. The forces driving their motion are determined by an internal (Le., intramolecular) harmonic potential obtained from the molecular dipole (linear) polarizability of water (1.444 A3) plus the external contributions. The intermolecular force on a charge follows from the electrostatic potential at the exact site where it is located in the solvent. These electrostatic potentials are calculated by using the familiar Ewald summation technique for charges under periodic boundary conditions, without any further modification.22 Once the fluctuating charges have been assigned a “mass parameter” m,, the equations of motion can be integrated by a standard molecular dynamics algorithm (in our case Verlet22). The mas mqhas no physical meaning and has the same function as the fictitious masses introduced in the constant-pressure and constant-temperature molecular dynamics methods.22 Although the values of these masses are in principle arbitrary, in practice they are chosen for technical convenience. The value we use here is m4 = where mH is the mass of a proton and a. is the atomic unit of length. This choice ensures that the polarization responds to the solvent motion on the time scale of the highest frequencies in the solvent, Le., of the librations, while the equations of motion can still be integrated with the standard time step of 2 . 5 fs. A well-known complication of the extension of charge site models with polarizable centers is the instability that can arise from a close approach of a positive point charge (proton) and an induced dipole moment. Some form of short-range damping of the induction is therefore necessary.I5 The main motivation for the construction of four point charges used here is that this damping can be achieved by simply broadening the fluctuating point charges into Gaussian charge distributions, mimicking in a very approximate way the finite size of the electronic density in the real molecule. The charge distribution at a tetrahedral site with position vector ri is thus written as
where qi is the total charge at the site i which fluctuates in time as outlined above. The short-range interaction in the TIP4P model is of the Lennard-Jones type A C u(r) = - - rI2 r6 with oxygen atoms as origin. For simplicity we adopt the TIP4P repulsive term. This leaves two adjustable parameters, $, of eq (22) Allen, M . P.; Tildesley, D . J . Computer Simulation of Liquids: Clarendon: Oxford, 1987.
Solvent Polarization and Hydration of CI-
The Journal of Physical Chemistry. Vol. 94, No. 16. 1990 6485
TABLE I: Parameters for Water-Water and Water-Ion Potentials QH, e QM, e rOM*
A
D E? A P,
10dAm, kJ AI2 m o P 10-3Cm, kJ A6 mo1-I 10-6ACI0,kJ AI2 mol-' 10-3Cclo, kJ A6 mol-I
TIP4P 0.52 -1.04 0.15 2.18
WKb 0.62 -1.24 0.15 2.60
2.510 2.552 16.53 6.114
5.000 4.850 9.328 2.999
PO1
0.4428 4 X -0.2214' 0.15 1.85 0.655 2.510 1.987 16.53 5.393
'See refs 3 and 7. bSee ref 5. ' Q M has been resolved in four fluctuating charges. The average (permanent) charge per site is quoted (see text and ref 14). dWidth of Gaussian distribution for fluctuating charges (see eq I ) .
1 and C of eq 2, which are specified by the following procedure.
First, we exploit the observation that the stability of the broad second maximum of the radial distribution function, goo(r), is very sensitive to the short-range damping parameter, f . This maximum is one of the most characteristic features of hydrogen bonding in water at ambient conditions and is only reproduced in the simulation when the value of f lies in a narrow range between 0.6 and 0.7 A. In the next step, the model is refined by using the experimental data on the total energy of the liquid.) In our dynamical scheme for the determination of the induced polarization, the comparison of the total energy to experiment requires some care. The total internal energy is the sum of the electrostatic interaction energy Ueland the polarization energy UP1 ( U ) = (Ud) + (Upd (3a) where the brackets denote a thermal average. U,, is the Coulomb interaction energy of all the charges in the system, including the fluctuating charges representing the polarization. For simplicity, in the following equation the finite width of these charge distributions is ignored. Thus, schematically
If we write UP, in terms of the induced dipole moments pi, we have
u,,.
=
c,Pi2 i La
The expression of pi in terms of the fluctuating charges can be found in ref 14. As was pointed out in ref 17 (see also refs 4, 5, and 14), the thermal average for UP, cannot be calculated by simply substituting the average induced dipole moments (pi) into eq 3c. The fluctuations (p: - (pi)2)(and in principle also the correlations) must be included to obtain the correct UPl. In our approach, the fluctuations are classical, and the polarization is in thermal equilibrium with the configurational degrees of freedom. This leads to a correction term of approximately )12ksT. In the actual simulation, of course, this problem does not arise because the polarization energy is evaluated as a time average of the complete expression of eq 3c for U . The point to note is that for our model the 3 / 2 k B Tterm must a d be added to the internal energy of the gas phase, even though the (pi) vanish for free molecules. With this caveat in mind, the simulation result for ( U ) can be compared to the experimental energy of evaporation. The value for C in the 0-0 potential of eq 2 is then found by adjusting it, with the previously specified ( held fixed, until the evaporation energy is in reasonable agreement with experiment (see Table 11). The final parameters used in the present study are listed in Table I . These values are similar to, but not identical with, the set of parameters of the first version of the model in ref 14. B. An Extended Mean-Field Model. The development of the extended mean-field model of ref 5 was motivated by consideration of several minor and major deficiencies of earlier mean-field
model~.l-~For most of these potentials the position of the first peak in the radial distribution function, g m ( r ) ,is not in agreement with the most recent improved neutron scattering data.*) More seriously, the dielectric constant of TIP4P (see ref 5) is underestimated, which may affect the ability of the model to solvate ions. The discrepancy for the dielectric properties can be attributed in part to the inaccurate representation of the molecular charge distribution and the fact that, with the exception of ref 4, the polarization energy is omitted from the expression for the internal energy. This inconsistency results in a reduced effective dipole moment and, consequently, a low dielectric constant. The mean-field potential of ref 5 has been designed with special attention to these three issues. The first step of the optimization procedure for the model started again from the experimental geometry. As with the TIP4P model,) a negative charge M, compensating the charges of the two protons, is placed along the molecular symmetry axis. The value of the proton charge QH and the separation roM between the oxygen (0)site and M is adapted to the experimental gasphase quadrupole tensor. The values of these parameters are given in Table I. The effective dipole moment of this array of charges turned out to be 2.6 D, in good agreement with the results of semianalytical mean-field calculations17 and information concerning polarization effects available in other s o ~ r c e s . ' ~ , ~ ~ A Lennard-Jones interaction between oxygen sites is also used in the extended mean-field model to describe the short-range interactions. In this case, both the coefficients A and C of eq 2 are treated as adjustable parameters. First, the range of possible values for A was narrowed down to a small interval by requiring that the position of the first peak in goo(r) coincides with the experimental value of ref 23. Then, the model was refined by an iteration process involving the pair of parameters A and C using agreement with the full experimental goo(r) and the internal energy as a guide. A further requirement was that the calculated pressure be small at the chosen density of 0.997 g cm-3 and temperature 25 O C . Special care was taken to account for the polarization selfenergy contribution to the total internal energy (see eq 3). For a mean-field model, this correction is no longer straightforward, since, unlike for the polarizable model, the polarization term of eq 3c is not available as an observable. Moreover, part of the polarization energy is contained in the electrostatic term of eq 3b, because it is precisely the choice of the enhanced charges that now determines the average induced dipole moment. The various contributions to the energy were separated following the scheme outlined in ref 17. The resulting parameters are given in Table I. Once again, the reader is referred to the original paper for more detail. C . Results f o r Pure Water. We have already reported preliminary results for the polarizable model,14 and the properties of the mean-field model were extensively reviewed in ref 5. For the present work we subjected the slightly modified polarizable model to further testing. The new results will be discussed in this section and compared to data taken from ref 5. As in the earlier ~ o r k , the ~ Jsample ~ studied consists of N = 216 water molecules in a box with periodic boundary conditions. For the polarizable model, the molecular dynamics trajectories were evaluated under (N,E,V) and (N,H,P) conditions. The runs in ref 5 were all (N,V,T) simulations. Figure 1 shows the results for g&), the oxygen-xygen radial distribution function. Although the polarizable model (pol) is able to reproduce the broad second peak with almost the correct maximum value, the maximum of the first peak is, however, too high. It is also shifted inward, a likely indication of too much attraction. It should be recalled that for the polarizable model the short-range repulsion parameter, A, was not optimized. For (23) Soper, A. K.; Silver, R. N. Phys. Reo. Len. 1982, 49, 471. Soper, A. K. Chem. Phys. 1984,88, 187. Soper, A. K.; Philips, M. G. Chem. Phys.
1986, 107, 47. (24) Coulson, C. A.; Eisenberg, D. Proc. R. SOC.London, A 1966, 291, 445. Barnes, P.; Finney, J. L. Faraday Discuss. Chem. SOC.1980, 69.
Sprik et ai.
6486 The Journal of Physical Chemistry, Vol. 94, No. 16, 1990
4, the induction is incorporated into the effective dipole moment of an array of discrete charges separated by a distance of order 1 A, a fact which likely enhances hydrogen bonding for a given dipole moment. This distinction also explains part of the rather large difference in polarization energy for the two models, which is almost a factor of 2 (see Table 11). Another effect that adds to Up, for the polarizable model is the finite temperature of the dynamical degrees of freedom associated with the polarization (see section IIA). As a consequence of fluctuations, the total induced dipole moment can also have an instantaneous component perpendicular to the permanent moment. A few relevant dynamical properties have been investigated in order to characterize the time scale of the relaxation processes in the liquid. Orientational relaxation was probed by means of the single-molecule time correlation function
4r
n
t.
0 0
CD
C(t)= (cos (ea(t).ea(0)))
0
2
4
6
B
l/A Figure 1. Oxygen-oxygen radial distribution function g,&) for water under ambient conditions. The dashed curves are experimental, taken from ref 23. The bold curves are calculated respectively for the extended mean-field model (WK) of ref 5 and the present fully polarizable model (pol). TABLE 11: Calculated Properties for Pure Water TIP4P WK POI expt 18.07 18.07 18.07 0.0 0.1 0.0 18.13 f 0.05 18.07 295 298 298 2.80 f 0.01 2.6d 2.18 2.60 24.9 13.2 -41.6 -42.5 -40.6 -41.5’ 3.3 f 0.5 1.1 f 0.3 1.5 f 0.4 2.48 3.6 f 0.2 8.8 f 0.3 10.8 f 0.5 3.6 f 0.3 10.1 f 0.4 7.9 f 0.4 2.7 f 0.3 6.8 f 0.3 7.8 0.4
*
Average internal pressure for isochoric system with volume given above. bAverage molar volume for isobaric system with pressure given above. eCompOnent of total dipole moment along molecular z axis. dSee ref 24. CCorrected for self-polarization energy; see section IIB and ref 5. /See ref 3. G e e ref 25.
the extended mean-field model (labeled WK), g m ( r ) is in far better agreement with experiment; here, also, the first peak is too sharp. Thermodynamic data are given in Table 11. To test the mechanical stability of the polarizable model a constant-pressure run was performed at an external pressure equal to the average virial pressure found from the constant-volume trajectory. The time average of the fluctuating volume agrees within statistical error with the constrained volume of the other run. One of the most interesting quantities that can be determined by means of the polarizable model is the effective dipole moment. Table I1 gives the average sum of the permanent and induced moment projected on the symmetry axis of each molecule. The quoted value of 2.80 D, which is significantly larger than the 2.6 D used in the extended mean-field (WK) mode, may seem somewhat surprising. However, comparable or even larger dipole moments have been reported recently for alternative polarizable model^.'^^^^ While the precise details of these alternative polarizable models are rather different, a common feature is that they are all based on a combination of distributed point charges and point polarizabilities. Thus, for this type of model, the condition of a stable hydrogen-bonded network seems to imply an induced (point) dipole of this magnitude (=I D). On the other hand, in the extended mean-field model, such as is used here and in ref
(4)
where ea is the unit vector which points along the a axis in the body-fixed coordinate frame of the molecule.6 In the present case, the z axis is the direction of the permanent dipole moment, the y axis is the vector connecting the two H atoms, and the x axis is perpendicular to the molecular symmetry plane. The characteristic times 7: determined from the slopes of In q ( t ) are given in Table 11. The results for the two models studied here are very similar. In a more detailed analysis in ref 5 it was concluded that the 7: are about a factor of 2 longer than those inferred, indirectly, from experimental data. The sluggishness of the long-time solvent dynamics is confirmed by the self-diffusion coefficients, also given in Table 11, along with experimental value.zs As is well-known, the TIP4P model3 and also the SPC modelZyield diffusion rates higher than experiment. In ref 5 it was argued that the increased rigidity can be understood as a consequence of the considerably stronger electrostatic forces in models that take the polarization energy into account. The same effect was noted in other recently proposed polarizable m0de1s.I~ 111. Hydration of the CI- Ion A . Ion- Water Interaction. In order to use the polarizable water model as a solvent, we need to adopt an ion-water interaction model. The ionsolvent potential is adapted from a set of models of which TIP4P is the member for pure water.’ The C1- ion is represented as a unit negative point charge with electrostatic coupling to the solvent. This interaction is supplemented with a Lennard-Jones term acting on the CI-O distance. We employ the same repulsive parameter, A (see eq 2), as in ref 7. The dispersion coefficient, C,is basically the same, except for a small reduction in order to adjust the ion-water parameter, C,to be consistent with the C coefficient for pure water. As explained in ref 14, the value of C for our polarizable model is smaller than the original TIP4P value. The ion-solvent potential used in combination with the WK model has the same general form. In this case, however, the A and C parameters have been optimized following a procedure similar to one for pure water outlined in section IIB and ref 5 . Both the location and minimum of the potential were modified to conform with information on the C1-0 distance and binding energy of the complex C1-.-H20.7 The final A and C values for the two models are listed in Table I. B. Results for the Anion Hydration. The molecular dynamics simulations of the CI- solvation were all performed with a system of 215 water molecules and a single anion. The time averages were determined over trajectories with a length varying between 30 and 40 ps. The technical implementation of the molecular dynamics code is the same as applied for pure water in the previous section. The oxygen-ion radial distribution function gc&) and the are shown in Figure corresponding function for hydrogen gCIH(r) (25) Krynicki, K.;Green, C. D.; Sawyer, D. W.Faraday Discuss. Chem.
Soc. 1918, 66, 199.
( 2 6 ) Mills,
R. G.Rev. Pure Appl. Chem.
1961, 1 1 , 79.
The Journal of Physical Chemistry, Vol. 94, No. 16, 1990 6487
Solvent Polarization and Hydration of CI-
3
Ci-pol. 2 -
-. n
3
b
E2
0 0
1
4
c
3
ci-np
A
2 1
0 0
2
0
6
4
Figure 2. Atom-atom radial distribution functions gcw(r) for the chloride anion solvated in water under ambient conditions calculated for the three different solvent models described in the text. The dashed and bold curves refer to X-H and X-0,respectively.
TABLE III: Calculated Properties for Cl- in Water TIP4P 7.2 f 0.1 7.0 f 0.1 3.27 f :Do? A 0.02 2.34 f r8IL A 0.02 3.60 f r8L A 0.04 A de8 2.70 f 0.05 A' D 2.18 &,, kJ mol-' -6 IO f 15 AHMl,kJ mol-' -370 f 20 AVMl, A' 43 f 21c D, IO4 m2 s-I 2.8 f 0.4 Tr, PS 7.5 f 1 no
WK
expt
POI
6.4 0.1 6.1 f 0.1 6.2 f 0.1 5.9 f 0.1 =6" 3.24 f 0.02 3.18 f 0.03 2.31 f 0.02 2.25 f 0.02 2.25 f O.0jb 3.70 f 0.04 3.50 f 0.05 3.2 f 0.05
3.0 f 0.06
5 f 5b
2.60 -670 f 15 -380 & 20
3.00 f 0.04 -539 f 2 -310 f 15 -30 f 40 1.1 f 0.3 13f2
-34W 39' 2.od =5
0.9 f 0.2 1612
2
4
6
rh
r h
"Estimate based on extrapolating the data.of ref 12 for ZnC12 to infinite dilution. bFrom the data of ref 12 for a 0.245 mol kg-l ZnCI, solution. CSeeref 7 and references therein. dFrom the data of ref 26 for NaCI extrapolated to infinite dilution.
2. The familiar characteristics of anion solvation are reproduced by all three models. Water molecules in the first solvation sheath are arranged in bond order with one of the two protons pointing toward the ion. The main structural difference between the three models is the number of molecules involved. The quantities no and nH in Table 111 are defined respectively by using the first minima of gclo and As suggested in the Introduction, the enhanced effective moment of the WK model, indeed, leads to the expulsion of a molecule from the first solvation shell. When the polarization is allowed to relax, the average dipole moment of molecules surrounding the ion increases by only a modest amount (less than IO%, see Table III), and a corresponding further small reduction of the coordination is achieved. The final hydration number for the polarizable model is less than that of the TIP4P model, an effect that can be attributed, at least in part, to an enhanced dipole moment. Table I11 also gives a few geometrical parameters: r position of the first maximum of gclx ( X is 0,H), and position of the second peak in gCIH, which arises from pendant H atoms of the coordinated molecules. The angle # is a measure of the deviation from perfect bond order. When = 0 the OH
k ik
+
Figure 3. Bold curves are neutron-weighted radial distribution functions AGcl(r) for the chloride anion solvated in water under ambient conditions calculated for the three different solvent models described in the text.
The dashed curves are experimental, taken from ref 12. bond, pointing from the ligand to the ion, is completely aligned with the vector, rWl, connecting the oxygen atom and the ion. The permanent dipole moment in this configuration is inclined by half the tetrahedral angle with respect to rWI. We verified that, as expected, the average induced dipole moment of the coordinated molecules is slightly rotated away from molecule symmetry axis in the direction of the bond vector, rWI. The calculated quantities of gclo and gClH can be compared to neutron scattering data. As is outlined in ref 12, the function AGcl obtained by the so-called "first-order difference" method is, in the limit of low concentration, a weighted average of g c o and gcw, namely (5) AGcdr) = &cio(r) 4- B g c d r ) Detailed expressions for the coefficients A and B, in terms of the appropriate coherent scattering lengths and concentration, can be found in ref 12. Figure 3 shows ACcl determined from radial distribution functions of Figure 2 along with the experimental datal2 in convenient units that compensate for solvent density effects. The first peak, in AGcl consists largely of the contribution from the first peak of &IH. For the polarizable model, the peak position is in reasonable agreement with experiment. The CIH distance for the mean-field models, TIP4P and WK, is somewhat larger (see also Table 111). The second feature in AGa is a sum of the first gclo peak and the second peak in gCIH, which nearly coincide. In the WK model these two contributions can be distinguished separately, whereas for the polarizable model gclo dominates. The energetics of the solvation is contained in Table 111. The absolute binding energy Ui,, of the ion with fluctuating solvent polarization is about 100 kJ mol-l less than for the mean-field models. Also, the total enthalpy of solvation AHsdv,which includes the solvent reconstruction energy, is smaller. For the polarizable model we also have made an attempt to evaluate the volume of solvation AV,, by comparing the average volume of a constant-pressure run to the average volume at the same pressure for the pure system (see Table 11). Although the accuracy is insufficient to draw definite conclusions, there seems to be a tendency to contraction-an observation that is contrary to the expansion observed in experiment. The last two entries in Table 111 give information on the dynamics of the solvated ion. The quoted diffusion coefficients confirm the picture already established for the pure liquid. (See the comments in ref 27.) While the diffusion for TIP4P is too
6488
The Journal of Physical Chemistry, Vol. 94, No. 16, 1990
fast for the WK extended mean-field model, the polarizable model diffusion is too slow. Note, however, that despite the error in absolute value, the diffusion coefficient of the ion is consistently about 20% smaller than the self-diffusion coefficient of the corresponding pure solvent (recall Table I I ) , a finding that is in excellent accord with experiment.26 Finally, in order to estimate the residence time for solvation sheath molecules, we have employed the definition of ref 6. There, the residence time T~ is obtained as the characteristic time for the decay of the correlation function N
P ( t ) = (CP,(t,O;t*))
(6)
’I I
Pj(t,t’;t*) is a property of molecule j which takes the value 1 if molecule j lies within the first coordination shell at both times t and t’and in the interim does not leave the coordination shell for any continuous period longer than t*. When these conditions are not satisfied, Pj(t,t’;t*) is zero. A molecule is counted as part of the first solvation shell when it is within a certain radius r, from the ion. In the present work, the parameter t* was set to 2 ps and rc was given the value corresponding to the position of the first minimum of gclo that was used to define the hydration number, no The results are listed in Table 111. The experimental residence time is generally estimated to be 5 ps or even less.s Hence, the T , for three models, including TIP4P, are far too long. This dynamical behavior is consistent with the enhanced structure in the radial distribution functions, AGcl (see Figure 3).
IV. Discussion and Conclusion Both the extended mean-field model and the polarizable model give a reduction of the coordination number of CI- in water, in comparison to the TIP4P potential. From this observation, we conclude that the repulsive interaction between the induced dipole moments of the molecules in the first solvation sheath is an important factor in determining the structure of the solvation shell. However, at least for the CI- ion, the enhanced effective dipole moment of the extended mean-field model is reasonably successful in accounting for these interactions. It remains to be determined whether or not this conclusion will be valid for other anions as well, such a F and Br- and the cations Na+ and Li+. A serious drawback of the type of models considered here is the apparent slow relaxation of the solvent configuration. This effect seems to be an inevitable consequence of including the polarization explicitly in the description of the solvent. Apparently, the increased electrostatic forces, which compensate the positive polarization self-energy, give the system too much rigidity. One may speculate that we are confronted here with the consequences of perturbing a fortuitous cancellation of errors in the original TIP4P model. The error, which we have attempted to remove, is the omission of the polarization self-energy. The other error, which still remains, is, perhaps, the approximation of a completely classical system of rigid molecules. Support for the importance of quantum librational motion can be found in the path integral studies of pure water,28where it was shown that the structural quantum effect amounts to an effective temperature increase of about 50 K. Clearly, a change of temperature of this magnitude would certainly accelerate the dynamics by the required amount! The consequence of using a purely empirical form for the pair potential should not be overlooked either. Further optimization of the potential will certainly be required, in particular realistic short-range atom-atom are interactions not given by A/rI2. Ideally, a polarizable model should be transferrable from the level of interaction in a dimer to the bulk without modification.16 For chloride ion-water clusters, a wealth of experimental data is (27) Friedman, H. L. Faraday Discuss. Cfiem.SOC.1988, 85, 1. (28) Kuharski, R. A.; Rossky, P. J. J . Cfiem. Pfiys. 1985, 82, 5164.
Sprik et al. available in literature.29 However, we have not yet subjected our model to this severe test. The reason is connected with our assumption of a finite temperature for the dynamical polarization fluctuations. This is not very realistic for a limited number of interacting molecules in a cluster. Also, the spurious finite temperature fluctuations of the polarization may act as a “friction” which contributes to the slow relaxation of the bulk fluid. The assumption of a finite temperature for the polarization is, however, not essential to our approach. In principle, it can be eliminated by “cooling” the system of fluctuating charges, while the temperature of the solvent configuration is kept constant at any desired value. With this additional dynamical constraint, the polarized model can, in fact, be considered as an application of the CarParrinello technique for electronic structure calculation^.^^^^' We are currently implementing this desirable extension of the simulation code. Other potentially serious simplifications in the modeling are the neglect of solvent intramolecular f l e ~ i b i l i t yand ~ ~ the polarization of the anion. It should be recalled that the polarizability of CI- is considerably larger than water. For an ion solvated in the bulk solvent, this approximation can be justified because, in the high-symmetry environment of the anion, the average induction on the ion is forced to vanish. For a cluster of irregular shape. the polarization of the anion must be included.33 Finally, for the case of two or more ions solvated in the bulk, the solute and solvent p ~ l a r i z a t i o nare ~ ~likely to be equally important. A fully self-consistent treatment of this example may perhaps be able to resolve the issue of the somewhat controversial predictions of mean-field models for the (Cl-)2 ion pair in water.27 Both simulations35and analytical integral equation calculation^^^ indicate that the potential of mean force exhibits an absolute minimum for the negative ions at contact, a configuration that is stabilized by the presence of bridging solvent molecules Hbonded to each anion.35 It is, perhaps, worth noting that the potential models employed in these calculation^^^ yield a hydration number of 10 for the CI- ion in dilute solution, a value that is clearly unreasonably large in view of the overwhelming contrary evidence presently available from experiment.I2 It will likely be instructive to obtain a quantitative estimate of the polarization contribution to the solvation characteristics of the (CI-)2 ion pair in water. Overall, the performance of the polarizable solvent model is, perhaps, a little disappointing, in view of its added computational expense. However, the main point that has emerged from the present study is that local polarization effects can alter hydration numbers of ions. This, in turn, will likely have ramifications in the modeling of more complex systems. Acknowledgment. We thank the referees for their careful reading of the manuscript. This research was supported by the National Institutes of Health (GM-40712) and the U S . National Science Foundation (CHE-872248 1 ). Computational resources were provided by the Stichting SURF from a grant out of “het National Fonds gebruik Supercomputers” in The Netherlands and the U S . National Science Foundation under CHE-8815130 at the University of Pennsylvania and DMR-88 19885 at the Pittsburgh Supercomputer Center. (29) Didzic, I.; Kebarle, P. J . Pfiys. Cfiem. 1970, 74, 1466. Arshadi, M.; Yamdagni, R.; Kebarle, P. J . Pfiys. Cfiem. 1970, 74, 1475. (30) Car, R.; Parrinello, M. Pfiys. Rev. Lptt. 1985, 55, 2471; 1988, 60, 204. (31) Sprik, M.; Klein, M. L. J . Cfiem. Pfiys. 1988, 89, 1592; 1989, 91, 5665. (32) Lie, G . C.; Clementi, E. Pfiys. Rev. A 1986, 33, 2679. Reimers, J. R.; Watts, R. 0. Cfiem.Pfiys. 1984, 91, 201. Toukan, K.; Rahman, A. Pfiys. Rev. B 1985, 31, 2643. ( 3 3 ) Foresman, .I.B.; Brooks, C. L. J . Cfiem. Pfiys. 1987, 87, 5892. (34) Goldman, S.; B a c k , P. Private communication, 1986. (35) Dang, L. X.; Pettitt, B. M. J . Am. Cfiem. SOC.1987, 109, 5531. J . Cfiem. Pfiys. 1987, 86, 6560. Preprint, 1989. (36) Pettitt, B. M.; Rossky, P. J. J . Cfiem. Pfiys. 1986, 84, 5836.
