1142
J. Phys. Chem. B 1997, 101, 1142-1147
Polarizable Ions in Polarizable Water: A Molecular Dynamics Study M. A. Carignano,† G. Karlstro1 m,‡ and P. Linse*,§ Physical Chemistry 1 and Theoretical Chemistry, Center for Chemistry and Chemical Engineering, Lund UniVersity, P.O. Box 124, S-221 00 Lund, Sweden ReceiVed: May 21, 1996; In Final Form: October 25, 1996X
The effect of the ionic polarizability on the solvation of positive and negative ions in water is investigated using molecular dynamic simulations. The NEMO many-body interaction potential for the water molecules is used, and a variety of values for the polarizability of the ions are considered. The increase of the polarizability leads to a larger electrical field at the ion through (i) a shrinking of the solvation shell around the ion and (ii) an increased probability of an asymmetric location of the ion in the cage. For the positive ions, whose polarizabilities are generally smaller than those of negative ions, these effects are visible only for the highest polarizability studied. However, for a given polarizability, the electrical field at an ion as well as its preference for an asymmetric location is larger for cations than for anions. The induced dipole moments in the first hydration layer are also investigated.
I. Introduction In spite of the long history of research on hydration of ions, the field is still very active and the theoretical understanding of such systems is certainly far from being complete. One reason for the difficulties in describing water and aqueous systems is the polarizability of the water molecule, which makes the many-body induction contribution to the forces and interaction energy important. In the homogeneous bulk phase, the effects of the polarizability are nearly isotropic and they can often be satisfactorily absorbed by a pairwise effective potential. However, when the symmetry is broken, like in the proximity of a surface or an ion, the explicit inclusion of many-body effects becomes important in order to obtain a realistic picture. Furthermore, in the study of hydration of ions, the interaction between the ion and the water molecules should also include the induction terms, depending on the polarizability of the ion. Only a few applications of simulation techniques to study hydration of ions with polarizable potentials are found in the literature, the reason being the requirements in computer time for calculating the induced dipoles, which were reduced to reasonable values with the computer technology of the late 1980s. A number of different ions at infinite dilution and ionwater clusters were studied with Monte Carlo1,2 and molecular dynamics.3-11 The most relevant work in relation to this paper is the MD simulation by Perera and Berkowitz,6 who have studied the effect of the sign of the charge and the polarizability of the ion in the cluster Cl--(H2O)20. In this paper we investigate the solvation of ions in water by using molecular dynamics simulations. In particular, we focus our attention on the effect of ionic polarizability on the solvation pattern. We consider both positive ions by looking at the solvation of Na+ and negative ions through the study of Cl-. However, it has to be stressed that the goal of this work is not to study the characteristics of the hydration of Na+ and Cl- as such but to understand the role of the ion polarizability in general. For that purpose, we have performed four different simulations for each ion with the only difference between them †
Physical Chemistry 1. Present address: Department of Chemical Engineering, Colburn Laboratory, University of Delaware, Newark, DE 19716. ‡ Theoretical Chemistry. § Physical Chemistry 1. X Abstract published in AdVance ACS Abstracts, January 15, 1997.
S1089-5647(96)01475-7 CCC: $14.00
being the ionic polarizabilities, which are varied around the normally adopted values for Na+ and Cl-, respectively. Normally smaller polarizabilities are observed for cations than anions: The former class varies from 0.028 Å3 for Li+ to 2.26 Å3 for Cs+, whereas the anion polarizabilities range from 1.31 Å3 for F- up to 7.4 Å3 for I-.12 From this fact, it is expected that the polarizabilities influence the local aqueous environment in the proximity of negative ions. This paper is organized as follows: The next section gives an outline of the intermolecular interactions and the methodology of the simulations as well as a description of the systems. Section III presents the results, and the conclusions are given in section IV. II. Potentials, Method, and Systems The many-body intermolecular potential used in these simulations for the interaction between the water molecules is the NEMO potential, recently developed and described in previous publications.13 The water molecule is assumed to have a rigid geometry with 104.52° and 0.9572 Å for the OHO angle and OH distance, respectively. Four point charges are present in the model, two of them are positive and they are located on the hydrogen atoms centers, whereas the other two are negative and placed in two virtual sites close to the oxygen. Detailed information on the water model is given in Table 1. This distribution of charge reproduces the dipole and quadrupole moments calculated at the ab initio SCF level.13 These point charges are allowed to interact with atom polarizabilities, which are also given in Table 1. The ion-water interaction is modeled as a Lennard-Jones pair potential plus the electrostatic and many-body induction terms. The total interaction energy is written as
U ) Uele + Uind + Urep + Udisp + ULJ
(1)
The first term is the electrostatic interaction
Uele )
∑i ∑j
qiqj rij
(2)
The summations run over all pairs of charges belonging to different molecules. The second term accounts for the many© 1997 American Chemical Society
Polarizable Ions in Polarizable Water
J. Phys. Chem. B, Vol. 101, No. 7, 1997 1143
TABLE 1: Water Molecule Geometry coordinate, Å
O
H1
H2
E1
E2
X Y Z
0.0 0.0 0.0
0.0 0.7570 0.5859
0.0 -0.7570 0.5859
0.1654 0.0000 0.2218
-0.1654 0.0000 0.2218
Charges and Polarizabilities q, e R, Å3
O
H
E
0 0.878 74
0.583 613 0.095 340
-0.583 613 0
Exchange Repulsion and Dispersion Parameters O-O O-H H-H
a, kcal/mol
b, Å-1
c, (Å kcal/mol)1/20
d, Å6 kcal/mol
99 774.35 5 220.48 40.92
3.954 3.992 2.033
2.494 1.297 0.000
156.22 46.61 8.64
body induction effects, due to the response of the polarizabilities to the electric field created by the rest of the system
Uind ) - 1/2
∑i BE0i ‚µbi
(3)
where the electrical field on atom i from the permanent charges is given by
B E0i
)
∑ j*i
qjb r ij (4)
rij3
and the induced dipole moment is
[ (
b µ i ) Ri B E0i +
3 ∑ j*i
)]
δij
rj b r ib r5ij
r3ij
b µj
(5)
Equation 5 together with eq 4 constitutes a set of coupled equations for the induced dipole moments. The third term in eq 1 is the repulsion between the water molecules
Urep )
∑i ∑j aije-b r
() cij
ij ij
+
20
rij
(6)
where the coefficients are given in Table 1. The fourth term of the total interaction energy is the dispersion energy between water molecules, which we write as
Udisp )
∑i ∑j
dij S(rij) r6ij
(7)
where the switching function S(r) is given by
S(r) ) 1 - e-(r/a)
n
(8)
with a being 1.2 Å and n ) 4. Finally, the remaining term in eq 1 is the Lennard-Jones interaction between the ion and the atoms of the water molecules
ULJ ) 4
ij ∑ i*j
[( ) ( ) ] σij
σij
12
-
rij
rij
6
(9)
All the Lennard-Jones parameters are given in Table 2. Molecular dynamics simulations in the canonical ensemble (NVT) were carried out using the MOLSIM14 package which includes the NEMO water-water potentials as one of the
Figure 1. Distribution function for the electrical field seen by the Na+ for different values of the ionic polarizability. The solid line corresponds to R ) 2.0 Å3, the dotted line to 0.4 Å3, the dashed line to 0.2 Å3, and the dot-dashed line to 0.1 Å3. The electric field’s unit is 106 statvolts/cm. (1 statvolt/cm ) 2.998 104 V/m.)
TABLE 2: Lennard-Jones Parameters , kcal/mol
σ, Å
0.02 0 0.12 0.02
3.3 0 3.8 3.2
O-Na H-Na O-Cl H-Cl
TABLE 3: Polarizabilities, Å3 real Na+ 12 simulated Na+ real Cl-1 12 simulated Cl-
0.148 0.1 3.76 2.0
0.2
0.4
2.0
3.7
5.5
7.3
installed potentials. The simulated system consisted of 216 water molecules and one ion in a cubic box with an edge of 18.6 Å. Periodic boundary conditions were adopted in the three Cartesian directions together with a spherical molecular cutoff of 8.5 Å. The classical equations of motion were integrated with a time step of 0.8 fs by the use of the velocity Verlet algorithm and a quaternion representation of the rigid molecules. The temperature was kept constant a 300 K by rescaling the velocities every time step. The induced dipole moments were calculated by solving eq 5 iteratively every fime time steps, and a first-order predictor was used to estimate their values in the intermediate steps. The validity of this procedure was addressed in ref 15. As was mentioned in the Introduction, we have performed four simulations for the Na+ ion and four simulations for the Cl- ion. The values of the polarizability of the ions in the different simulations are given in Table 3. They cover the ranges representative of alkali and halide ions. Again, it is only the polarizabilities which differ in the two sets of simulations. Any differences in the results within a series can thus be attributed to different polarizabilities. We started the simulations for each ion with the smallest value of the polarizability, running an equilibration of 40 ps and a production run of 50 ps. III. Results We begin this section by presenting the results obtained for the systems containing the positive ions. The distribution functions for the electrical field seen by the Na+ are shown in Figure 1 for all the different values of the polarizability, whereas the mean values of these distributions are given in Table 4. For the three smaller polarizabilities the distributions are almost identical, but for the largest one an important shifting toward larger values of E is observed. The electrical field at the ion is the result of an interplay among the charges and polarizabilities of the atoms in the water molecules and those of the ion itself. In a first approach, assuming that changes in the ionic
1144 J. Phys. Chem. B, Vol. 101, No. 7, 1997
Carignano et al.
TABLE 4: Results for the Cationic Systems polarizability, Å3
E, 106 statvolts/cm
rNa-O, Å
njONa
0.1 0.2 0.4 2.0
0.183(0.002) 0.185(0.002) 0.199(0.002) 0.342(0.004)
2.41 2.41 2.41 2.38
5.72(0.02) 5.71(0.02) 5.72(0.02) 5.48(0.02)
polarizability do not affect significatively the structure of the liquid, it can be concluded from eq 5 that an increase in the polarizability will be reflected in a larger electrical field at the ion. However, in order to have a better understanding of the problem, a more careful analysis considering the effects of the ionic polarizability on the structure of the liquid around the ion should be done. The ion-oxygen radial distribution function (rdf) gNa-O(r) and the corresponding ion-hydrogen rdf gNa-H(r) are shown in Figure 2 for the four different values of the polarizability considered for Na+. The first thing to note is that there is almost no difference between the rdf’s corresponding to the four different cases. Only for the largest polarizability is the first maximum of gNa-O(r) shifted slightly toward the ion (see Table 4). The pattern shown is the characteristic one for the hydration of a small cation: (i) a well-defined first solvation layer exists (a sharp peak in the ion-oxygen rdf followed by a deep minimum) and (ii) the water molecules have their oxygens pointing toward the ion, leaving their hydrogens to face the bulk water (the ion-hydrogen rdf has its peak at 0.5 Å further away from the maximum of the ion-oxygen rdf). The coordination number njONa, defined as the average number of water molecules within a distance up to the first minimum of the gNa-O(r), decreases slightly for the largest polarizability as is shown in Table 4. These results are in good agreement with experimental data obtained by X-ray diffraction,16 which give 2.4 Å for the distance Na-O and 5.0 for the coordination number. Additional and complementary information on the structure of the liquid around the ion can be obtained from the angular distribution function (adf) P(cos θ). We define θ as the angle formed between the vector connecting the ion and the water center of mass and the direction of the water dipole moment. In Figure 3 we show P(cos θ) for the water molecules in the first solvation shell, i.e., those molecules that have their centers of mass closer than 3.20 Å (3.17 Å for R ) 2.0 Å3) from the ion. As we can see, P(cos θ) has a maximum at θ ) -180° which corresponds to the orientation of the water σ lone pair toward the Na+. However the maximum is quite broad, in agreement with the width of the peak of the hydrogen-ion rdf, indicating large fluctuations around the equilibrium conformation. As was pointed out for the rdf’s, only small changes are observed in the adf for the different polarizabilities considered here. To complete our structural analysis of the infinitely dilute ionic solution, we study the relative position of the ion with respect to the surrounding water molecules. In order to do that, we define the quantity ri to measure the distance between the ion and the center of the cavity according to
ri ) |〈b r w〉 - b r ion| )
|∑
r w(j) jW(j)(b
|
-b r ion)
∑jW(j)
(10)
rion|-6, and The weighting factor W(j) was set equal to |r bw(j) - b rion| ) 7 Å. To understand a cutoff was imposed at |r bw(j) - b the meaning of this quantity we consider some particular cases. If the ion is located in the geometrical center of a spherical cavity in the water, 〈r bw〉 will be the position of the (weighted) center of mass of the water molecules within a sphere of radius
Figure 2. Ion-oxygen (gNa-O(r)) and ion-hydrogen (gNa-H(r)) radial distribution functions for the different cationic systems. Lines are as in Figure 1.
Figure 3. Angular distribution functions for water molecules in the first hydration shell for the different cationic systems. Lines are as in Figure 1.
Figure 4. Distribution function for ri, a measure of the asymmetric location of the ion in the solvation cage, for the different cationic systems. Lines are as in Figure 1.
7 Å around the ion. Therefore, 〈r bw〉 will be at the center of the cavity and then ri ≡ 0. On the other hand if the ion is in a different position, the average 〈r bw〉 will still be close to the center of the cavity, so that ri > 0. In Figure 4 the distribution function P(ri) is plotted for all the simulated cationic systems. We see that for the three lower polarizabilities there is no difference in the distributions. However, for the highest value of the polarizability the distribution is shifted toward larger values of ri. Considering that the main features of the rdf’s are the same when increasing the ion polarizability up to 2.0 Å3, the cage should be essentially the same for the four different cases. Then, the origin of this shifting in P(ri) has to be in a higher affinity of the ion for the periphery of the cage. In other words, the ion wants to avoid
Polarizable Ions in Polarizable Water
J. Phys. Chem. B, Vol. 101, No. 7, 1997 1145
Figure 5. Distribution function for the induced dipole moment of the water molecules in the first solvation shell for the different cationic systems. The lines with circles correspond to the idm of the water molecules in the bulk. The other lines are as in Figure 1.
Figure 7. Ion-oxygen (gCl-O(r)) and ion-hydrogen (gCl-H(r)) radial distribution functions for the different anionic systems. Lines are as in Figure 6.
TABLE 5: Results for the Anionic Systems
Figure 6. Distribution function for the electrical field seen by Cl- for different values of the ionic polarizability. The solid line corresponds to R ) 7.3 Å3, the dotted line to 5.5 Å3, the dashed line to 3.7 Å3, and the dot-dashed line to 2.0 Å3. The inset shows the results corresponding to R ) 2.0 Å3 for Na+ (open symbols) and Cl- (filled symbols).
the geometrical center of the cage. This effect is not clearly visible from the rdf (see Figure 2), but it is amplified in eq 10 by the strongly distance dependent weighting factor. This behavior can be understood from classical electromagnetism.17 If we consider a polarizable point charge in a cavity with radius R in a continuous dielectric medium with an infinite dielectric permittivity, the field seen by the polarizability is proportional to -qr/(R3 - R), r being the distance from the center of the charge. Therefore, the polarizable charge wants to move toward the dielectric. Since the energetic contribution originated from the coupling between the field and the polarizability is proportional to the field squared and the polarizability, the effect is larger for larger polarizabilities. The induced dipole moment (idm) of, or the electrical field on, the water molecules will also be influenced by the presence of the ion. Figure 5 shows the distribution functions for the induced dipole moment of the water molecules in the first solvation shell as well as that for the water molecules in the bulk. The idm is always smaller in the first solvation shell than in the bulk. The increase in the polarizability of the ion leads to a very small increment in the idm of the neighboring water. The idm of molecules in the second hydration layer (not shown in Figure 5) do not show any influence from the presence of the ion. We will now concentrate on the anionic systems and their solvation properties. As was mentioned before, negative ions have larger polarizability that positive ions, and we have chosen four representative values covering that range, including the real polarizability for Cl- (see Table 3). Figure 6 shows the distribution functions for the electric field on the site occupied
polarizability, Å3
E, 106 statvolts/cm
rCl-O, Å
O njCl
H njCl
2.0 3.7 5.5 7.3
0.173(0.002) 0.248(0.003) 0.435(0.004) 0.976(0.005)
3.42 3.41 3.36 3.21
8.44(0.03) 8.34(0.03) 7.86(0.03) 7.69(0.04)
7.08 6.79 5.91 5.11
by the ion from the four different simulations with Cl-. In contrast with the previous case, all the curves are distinctly different and show a shift toward larger values of E with increasing the polarizability. The mean values of the distribution are given in Table 5. A question that arises concerns the relation between the electric field on the Na+ and the Cl- for the same value of the polarizability. The inset of Figure 6 shows these distributions for R ) 2.0 Å3, and it shows that the electrical field on the cation is twice as large as that on the anion. The reason for this relation has to be found by looking at the differences between the systems. Na+ is smaller than Cl-, so that the water molecules are closer, in that case favoring a higher electrical field. Secondly, the cation has the oxygens of the water molecules as the first neighbors, whereas the anion has the hydrogens. The oxygen atom has a much larger polarizability and a larger charge than the hydrogen atom, so again this effect enhances the electrical field on Na+. Figure 7 shows the rdf’s gCl-O(r) and gCl-H(r) for the four different polarizabilities. First we note that the curves for different polarizabilities are not similar as in the case of Na+ (see Figure 2). Here the polarizability has a more pronounced effect. The position of the first maximum is gCl-O(r) is located at 3.42 Å for the smallest polarizability (2.0 Å3) and is reduced to 3.21 Å for the highest one (7.3 Å3). Correspondingly, the first maximum in gCl-H(r) shifts from 2.51 to 2.28 Å. This effect is reflected also in the coordination numbers, which are summarized together with the other coordination properties in Table 5. These numbers are larger than experimental results obtained by neutron diffraction technique;18 those results give 2.25 ( 0.05 Å for the Cl-D distance and a hydration number smaller than 6. The first solvation layer is not as well-defined for these anions as for the cations. A weak minimum can be distinguished for the three smaller polarizabilities, but no minimum at all is present for R ) 7.3 Å3. In order to calculate the idm and the adf for the water molecules in the first solvation shell, we choose to use the minimum in the gCl-O(r) for the
1146 J. Phys. Chem. B, Vol. 101, No. 7, 1997
Carignano et al.
Figure 8. Angular distribution functions for water molecules in the first hydration shell for the different anionic systems. Lines are as in Figure 6. Figure 10. Distribution function for the induced dipole moment of the water molecules in the first solvation shell of the different anionic systems. The lines with circles correspond to the idm of the water molecules in the bulk. The other lines are as in Figure 6.
Figure 9. Distribution function for ri, a measure of the asymmetric location of the ion in the solvation cage, for the different anionic systems. Lines are as in Figure 6.
three lower polarizabilities and for the highest one where no minimum occurs we substract the difference in the radius for the maximum of the gCl-O(r) for the two higher polarizabilities from the radius used for the second highest polarizability (see Table 5). The shrinking of the first hydration shell with increasing the polarizability was also observed by Perera and Berkowitz for a cluster of a Cl- ion and 20 water molecules.6 They carried out two simulations with R ) 0.0025 and 3.25 Å3, respectively. They found a change in the coordinate number from 5.5 for the small polarizability to 4.5 for the large one. Figure 8 shows the adf’s for the water molecules in the first solvation layer. The maximum in the distribution is centered around 0.61, corresponding to an angle θ ) 52°. This reflects a preference for the formation of linear hydrogen bonds with the Cl-, a preference that is slightly enhanced as the polarizability is increased from 2.0 to 7.3 Å3. The position of the anion in the cage also shows a more dramatic dependence on the polarizability as compared to the cation. Figure 9 shows the distribution P(ri) for the four simulations with Cl-. The curves show a gradual displacement of the maximum with increasing polarizability, showing a tendency of the ion to avoid the center of the cage as was observed already for the cations (Figure 4). However, since the anions have larger polarizabilities, this effect is now clearly evident. At this point we have to note the qualitative similarity between the distribution for the electrical field at the ion and that of the position of the ion in the cage, i.e., Figures 1 and 4 for the Na+ and Figures 6 and 9 for the Cl-. Our results show a clear correlation between both quantities. As was mentioned before, both properties should respond in the same direction. In Figure 10 we display the distribution functions for the idm of the water molecules in the first solvation layer for the different anionic systems together with those of the water in the bulk. Again as the polarizability is increased, now from 2.0 to 7.3 Å3, the distribution is shifted to larger idm and is broadened. For the smaller values of R, the idm of the first shell is an average lower than that in the bulk, similarly to the case of
Na+. However, at increasing polarizability a crossover is observed at R = 5.5 Å3. The broadening of the distribution for the highest R is consistent with the fact that the environment of the water molecules in the first hydration shell becomes more heterogeneous due to a more asymmetric anion location. Finally we note that the long-range interaction in ionic and dipolar systems has to be approximated in molecular simulations and that the results consequently are influenced by the method used. This has been a subject of long-standing interest, but still no ultimate solution is in sight.19,20 The two approaches mostly used are (i) the Ewald summation and (ii) the minimum image or the related spherical cutoff. Both methods have their weaknesses, and the question which of them to prefer in liquid systems, if any, is not settled yet. IV. Conclusions We have presented a detailed study of the relation between the structural properties of an infinitely dilute ionic solution and the polarizability of the ion using molecular dynamic simulations. We have used a polarizable model to describe the system, and we have considered a number of cases with different ionic polarizability both for anions and cations. The electrical field at the ion increases when increasing the polarizability. The most important effects are observed for the anions, due to their larger polarizabilities. The relation between the electric field and the polarizability is nonlinear, but the effect is irrelevant for the three smaller polarizabilities considered for the cations in this work. Also, a large value of the polarizability favors an asymmetric location of the ion in the solvation cage, and a correlation between electrical field at the ion Ei and the weighted position of the ion in the cage ri was found. The polarizability also affects the structure of the water molecules in the first solvation layer. The main effect is a shrinking of the solvation cage for high polarizability, as was already observed by Perera and Berkowitz in clusters of a Clion and 20 water molecules. The orientation of the water molecules in the first solvation layer does not show an important influence from the polarizability. The induced dipole moment of the water molecules in the first hydration layer also displays some sensitivity to the value of the polarizability of the ion, but this is relevant only for anions. For the same polarizability, all the effects are larger for the cations than for the anions. The reason is that the former are smaller in size, so that all the interaction with the closest water molecules are stronger. Another reason is the sign of the ionic charge, namely, for the cations the oxygens of the water molecules are first neighbors and they have a larger charge and polarizability than the hydrogens. The anionic polarizabilities
Polarizable Ions in Polarizable Water vary over a larger range and have values larger than the corresponding cationic polarizabilities. The present results suggest that the surface activity of ions increase with their polarizabilities. This is well in line with the lyotropic series.21 Thus, we believe that the main origin of the lyotropic series is the variation of the polarizabilities (partly in relation to the size) at the ion. Finally, specific simulations for Cl- and Na+ ion with new NEMO potentials for the ionwater interactions are currently underway. In this later study we will also consider the hydration of the ions in small clusters. Acknowledgment. M.A.C. acknowledges the Swedish National Science Research Council (NFR) for a postdoctoral fellowship. References and Notes (1) Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1985, 83, 2923. (2) Cieplak, P.; Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1987, 86, 6393. (3) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J. Am. Chem. Soc. 1991, 113, 2481. (4) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1991, 95, 1954.
J. Phys. Chem. B, Vol. 101, No. 7, 1997 1147 (5) Dang, L. X. J. Chem. Phys. 1992, 96, 6970. (6) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1992, 96, 8288. (7) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1993, 99, 4222. (8) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1994, 100, 3085. (9) Sremaniak, L. S.; Perera, L.; Berkowitz, M. L. J. Phys. Chem. 1996, 100, 1350. (10) Dang, L. X.; Garrett, B. C. J. Chem. Phys. 1993, 99, 2972. (11) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 101, 7873. (12) Pyper, N. C.; Pike, C. G.; Edwards, P. P. Mol. Phys. 1992, 76, 353. (13) Wallqvist, A.; Ahlstro¨m, P.; Karlstro¨m, G. J. Phys. Chem. 1990, 94, 1649. (14) Linse, P.; Wallqvist, A. MOLSIM 1.3; Lund University: Lund, Sweden, 1992. (15) Ahlstro¨m, P.; Wallqvist, A.; Engstro¨m, S.; Jo¨nsson, B. Mol. Phys. 1989, 68, 563. (16) Skipper, N. T.; Neilson, G. W. J. Phys.: Condens. Matter 1989, 1, 4147. (17) Bo¨ttcher, C. J. F. Theory of Electric Polarization; Elsevier: Amsterdam, 1973. (18) Powell, D. H.; Barnes, A. C.; Enderby, J. E.; Neilson, G. W.; Salmon, P. S. Faraday Discuss. Chem. Soc. 1988, 85, 137. (19) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, U.K., 1987. (20) NRCC Proc. No. 9; Ceperly, D., Ed.; National Technical Information Service: Springfield, VA, 1980. (21) Collins, K. D.; Washabaugh, M. W. ReV. Biophys. 1985, 18, 323.
