Chemical Physics Letters 393 (2004) 249–255 www.elsevier.com/locate/cplett
Accommodation coeﬃcients for water vapor at the air/water interface John Vieceli *, Martina Roeselov a, Douglas J. Tobias
Department of Chemistry, University of California, Irvine, CA 92697-2025, USA Received 18 May 2004; in ﬁnal form 11 June 2004 Available online 2 July 2004
Abstract The accommodation of water vapor at the air/water interface is studied using classical molecular dynamics computer simulations. From 250 initial conditions of a gas phase water molecule approaching the surface of liquid water with a thermal impact velocity, the thermal and mass accommodation coeﬃcients are determined to be 1.0 and 0.99, respectively, at 300 K. When the approaching above the liquid, it is accelerated towards the surface, which is followed by a slight rotational excitation upon water molecule is 5 A impact. Within 4 ps of striking the surface, the total kinetic energy is equilibrated and the equilibrium number of water molecules in the ﬁrst and second solvent shells is obtained. 2004 Elsevier B.V. All rights reserved.
1. Introduction The process of gas uptake at a liquid surface can be characterized by two probabilities of accommodation. The thermal accommodation coeﬃcient, S, is a measure of the number of gas molecules that equilibrate to the temperature of the liquid Number of molecules equilibrated to liquid temperature S¼ : Number of molecules impinging on liquid surface ð1Þ
The mass accommodation coeﬃcient, a, describes the probability of the gas molecule being incorporated into the liquid a¼
Number of molecules absorbed into liquid : Number of molecules impinging on liquid surface ð2Þ
An accurate determination of these probabilities is crucial to understanding the nucleation and growth kinetics of cloud droplets and, subsequently, to cloud albedo and climate change [1,2]. The accommodation process is treated commonly as a two-step mechanism * Corresponding authors. Fax: +1-949-824-8571 (J. Vieceli), +1-949-824-9920 (D.J. Tobias). E-mail addresses: [email protected]
(J. Vieceli), [email protected]
0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.06.038
consisting of gas adsorption at the liquid surface followed by desorption into the gas phase or solvation in the bulk liquid . The rate constant for gas adsorption at the liquid surface is proportional to S, while a is related to the rate constants for desorption into the gas phase, kdesorb , and solvation in the bulk liquid, ksolv , by a¼
ksolv : ksolv þ kdesorb
Unfortunately, various experimental determinations of S and a for water vapor on liquid water are inconsistent. Reported values of S range from 0.7  to 1.0 [5,6]. The range of experimentally determined a values is even larger than that for S. A review article on this subject shows that a values in the range from 0.03 to 1.0 have been reported . Subsequent to that review article, the temperature dependence of a has been studied using a droplet train apparatus . A negative temperature dependence is observed with a ¼ 0:32 at 258 K and a ¼ 0:17 at 280 K. Molecular dynamics (MD) computer simulations provide another method for the determination of accommodation coeﬃcients. This approach has been used previously for the determination of a, but a determination of S from MD simulations is not available to our knowledge. A recent study by Morita et al.  reports a > 0:99 at 273 K. While the MD simulations are most consistent with the experimental
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255
conditions in the droplet train apparatus, the values of a obtained from these methods are not in agreement. Another MD simulation that focused on the temperature dependence of a obtained values in the range from 0.961 (T ¼ 330 K) to 0.286 (T ¼ 550 K) . It is evident from the lack of agreement between experiments and between experiments and simulations that further work is necessary on the accommodation of water vapor on liquid water. Reconciliation of experimental and simulation results is especially compelling in light of the importance of accurate S and a values to atmospheric chemistry [1,2]. In this Letter, results from MD simulations for the thermal and mass accommodation coeﬃcients of water vapor on liquid water at 300 K are presented. The simulation temperature is selected to be intermediate to previous MD simulations at 273 K  and 330 K . The objectives of this study are to develop an understanding of the accommodation process from a molecular perspective on a picosecond time scale and to determine values for S and a at 300 K. S and a are calculated according to Eqs. (1) and (2), respectively, by preparing 250 initial conditions of a gas phase water molecule with a thermal velocity approaching the surface of liquid water and monitoring the incoming molecule’s velocity and position. In addition to water [8,9], this approach has been used previously to study the mass accommodation of ethanol , methanol  and HO2  on liquid water and hydroxyl and ozone on aqueous salt solution surfaces . Using an alternative approach, Taylor and Garrett  calculate the potential of mean force for the transfer of a solute from the gas phase to the bulk liquid and use transition state theory to determine a according to Eq. (3) for ethanol and ethylene glycol on liquid water. Both approaches yield a consistent result for the mass accommodation of ethanol at the air/liquid interface of water and should be used in conjunction so that the short and long time aspects of the accommodation process are studied . From the simulation results presented in this Letter, it is determined that S ¼ 1:0 and a ¼ 0:99 at 300 K. Although a couple of water molecules directly scatter upon impact with the surface, the water molecules are equilibrated to the liquid temperature. As the incoming molecule approaches, the translational kinetic energy in the direction normal to the surface increases. Upon striking the surface, the rotational kinetic energy is slightly excited and the total excess kinetic energy is dissipated within 4 ps. The structure in the ﬁrst and second solvent shells of the adsorbed water molecule is equilibrated within 4 ps of striking the surface also. Desorption from the surface is not observed within 90 ps, which is consistent with the potential of mean force for the transfer of a water molecule across the air/liquid interface of water [14,15].
The remainder of this Letter is organized as follows. In Section 2, the MD methodology is described. The results and discussion are presented in Section 3, followed by conclusions in Section 4.
2. MD methodology Conﬁgurations of the liquid to be used for the accommodation study are generated from a simulation of 864 water molecules in a slab geometry with box di (x) · 30 A (y) · 100 A (z). The z dimensions of 30 A mension is perpendicular to the air/liquid interface. Following a 500 ps equilibration, a 250 ps equilibrium simulation is performed. Conﬁgurations from the equilibrium simulation, stored every 50 ps, are used to generate ﬁve initial conﬁgurations for the accommodation study. For the remainder of this Letter, the term air is used in reference to the vacuum region of the simulation cell. It has been shown previously that this simulation cell is a reasonable model for a small patch of the air/water interface based on the density of nitrogen in air and the free volume of the simulation cell . For the accommodation study, a single gas phase above the water molecule is placed approximately 15 A The veair/liquid interface in the xy-plane at z ¼ 30 A. locity of each atom is assigned from a Maxwell–Boltzmann distribution at 300 K, subject to the constraint that the component of the center of mass velocity normal to the air/liquid interface is directed towards the surface. This procedure is repeated ﬁfty times for each of the conﬁgurations from the neat simulations by placing the gas phase water molecule at diﬀerent grid points in the xy-plane, which generates a total of 250 trajectories. The trajectories are run for a total of 90 ps. During the ﬁrst 15 ps, the velocity and position of the incoming water molecule are stored every 0.1 ps for a detailed kinetic energy analysis that can be correlated with the position. During the ﬁnal 75 ps, the position is stored every 1.0 ps. The water molecules are described using the three-site water model of Caldwell and Kollman  that explicitly includes the induced polarization (POL3). The enthalpy of vaporization, density, diﬀusion constant and bulk liquid structure that are calculated using this model exhibit good agreement with experimental data. The POL3 water model also provides a reasonable description of the surface of liquid water. For example, the surface tension of POL3 water at 300 K (55 ± 6 mN/m) is within 23% of the experimental value at 303 K (71 mN/m) . (A diﬀerent polarizable water model yields a surface tension value of 92 ± 5 mN/m at 298 K .) For POL3 water molecules located in the surface region between 90% and 10% of the liquid water density, the percentage of OH bonds that are not involved in hydrogen bonds is 27%. This agrees with a lower bound of
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255
For the process of a thermal gas phase molecule impinging on a liquid surface, there are several possible outcomes for the incident molecule. The molecule may directly scatter back into the gas phase or adsorb at the interface. If adsorbed, the molecule, at a later time, may desorb back into the gas phase, remain adsorbed at the interface or be absorbed into the bulk liquid. The number of trajectories that result in each of these outcomes is determined from the data shown in Fig. 1. In Fig. 1a, the oxygen atom z position of the incident water molecule is shown for all 250 trajectories as a function of time. The air/water interface, which is deﬁned by the region between 10% and 90% of the liquid water density for this Letter, is located between the dotted horizontal lines. From the 250 initial conditions, 248 of the trajectories adsorbed at the interface and either remained at the interface or dissolved into the bulk liquid. Desorption from the surface following adsorption is not observed on the 90 ps simulation time scale. Two of the trajectories result in direct scattering events upon impact with the liquid surface. In Fig. 1b, the probability distributions of the oxygen atom z position for the 248 adsorbed trajectories at 30, 60 and 90 ps are shown. The population in the bulk region increases as time progresses due to diﬀusion of water molecules from the surface region. The percentage of these trajectories that are located in the bulk region is 66%, 75% and 84% at 30, 60 and 90 ps, respectively. The predicted outcome of the 16% that are located in the interface at 90 ps, which is relevant to the value of a, is
3. Results and discussion 3.1. The overall process
20% estimated from the vibrational sum frequency generation spectrum reported by Du et al. . The induced dipoles at each time step are determined self-consistently using the standard iterative procedure [21,22]. In order to increase computational eﬃciency, the internal geometry of each water molecule is held rigid in the equilibrium conﬁguration using the SHAKE algorithm . The omission of molecular vibrations from the simulations is expected to have a negligible eﬀect on the results since vibrational excitations are unlikely to be caused by thermal collisions. All of the simulations are performed in the NVE ensemble with an integration time step of 1.0 fs at 300 K. Three-dimensional periodic boundary conditions are used. Electrostatic interactions are calculated according to the smooth particle mesh Ewald method . The real space part of the Ewald sum and the van der Waals interac The simulations are pertions are truncated at 12 A. formed using version 6 of the AMBER software package .
0 -10 (b)
10 z (Å)
Fig. 1. (a) The oxygen atom z position of the incident water molecules as a function of time for those that adsorbed (solid lines) and directly scattered (dashed lines). The dotted horizontal lines and labels to the right of the panel denote the regions of the system, as explained in the text. (b) Probability distributions of the oxygen atom z position for the 248 adsorbed trajectories at 30 (dashed line), 60 (thin solid line) and 90 (thick solid line) ps. The dotted vertical lines and labels denote the regions of the system.
discussed in Section 3.3. Prior to addressing the value of a from these simulations, the process of thermal accommodation is analyzed and the value of S is determined. 3.2. Thermal accommodation The thermal accommodation process is studied by calculating a time-dependent function that determines whether each incident water molecule reaches thermal equilibrium with the liquid. This function is given by dKEtot ðtÞ ¼ KEtot ðtÞ hKEtot i;
where KEtot ðtÞ is the time-dependent total kinetic energy of the incident water molecule and hKEtot i is the equilibrium ensemble average of the total kinetic energy for a water molecule with three internal constraints at 300 K (hKEtot i ¼ 3kB T , where kB is Boltzmann’s constant). Based on the previous study by Wilson and Pohorille  on the accommodation of ethanol, this function is expected to rise above zero as the gas phase molecule approaches the surface, followed by a relaxation on the
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255 3
density (molecules/Å ) 0
4 tphase (ps)
Fig. 3. The oxygen atom z position of the incident water molecules as a function of time for the average of those that adsorbed (thick solid line) and for those that directly scattered (dashed lines). The water density proﬁle in the direction normal to the air/liquid interface (thin solid line) is also shown. The lower x-axis refers to the time passed and the upper since the incident water molecule reached below 21 A x-axis corresponds to the water density.
4 tphase (ps)
picosecond time scale, and then oscillate around zero at thermal equilibrium. Since the magnitude of the velocity component normal to the interface is randomly sampled for each trajectory, the gas phase water molecules reach the surface at diﬀerent times. In order to average or compare the time-dependent results from diﬀerent trajectories, the time for all of the trajectories must be phased to a common point. This point is selected to be the time step when the gas phase molecule ﬁrst passes below z ¼ 21 A, which is approximately 5 A above the liquid surface. The time axes in the remaining ﬁgures to be presented refer to the time passed since crossing the xy plane at and are denoted by tphase . z ¼ 21 A In Fig. 2a, the average of dKEtot ðtÞ for the 248 adsorbed trajectories and dKEtot ðtÞ for the two directly scattered trajectories is shown. This data can be correlated with the position of the incident molecules using the data in Fig. 3. The oxygen atom position as a function of time for the average of the adsorbed tra-
Fig. 2. (a) The function deﬁned by Eq. (4), dKEtot ðtÞ, for water molecules that adsorbed (solid line) and directly scattered (dashed lines). The time axis refers to the time passed since the approaching water (b) The distributions of molecular speed molecule reached below 21 A. at tphase ¼ 3 (dashed line), 4 (thin solid line) and 5 (thick solid line) ps for the entire ensemble of trajectories and the Maxwell–Boltzmann distribution of molecular speed for water at 300 K (dotted line).
jectories and for the two directly scattered trajectories is shown in Fig. 3, in addition to the water density proﬁle in the z dimension. For the adsorbed trajectories, the total kinetic energy rises due to an attractive potential energy with the liquid surface and reaches a maximum at 0.6 kcal/mol above the equilibrium value. This corresponds to a heating from 300 K to approximately 400 K in one picosecond. From the average position of the adsorbed trajectories in Fig. 3, the peak in dKEtot ðtÞ at 1.0 ps is correlated with the incoming molecule striking the liquid surface. The excess kinetic energy is dissipated within 2 ps of striking the surface and then oscillates around zero at thermal equilibrium. During these oscillations, the molecules are located on average at a position that corresponds to 90% of the bulk liquid density. The amount of excess kinetic energy and the time scale for energy dissipation are similar to previously reported results for ethanol accommodation on the surface of liquid water at 310 K . The behavior of the two trajectories that resulted in direct scattering is quite diﬀerent from the average for the adsorbed trajectories. In particular, both molecules undergo a more highly energetic collision with the liquid surface than the adsorbed trajectories and, then, directly scatter back into the gas phase.
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255
0.8 x translation
In order to assess the thermal accommodation coefﬁcient, the distributions of molecular speed at tphase ¼ 3, 4 and 5 ps for the entire ensemble of trajectories are compared to the Maxwell–Boltzmann (MB) distribution of speed for water at 300 K in Fig. 2b. While the distributions at tphase ¼ 3 and 4 ps exhibit slight deviations from the MB distribution near the peak, the distribution at tphase ¼ 5 ps agrees well with the MB distribution. This indicates that the incident water molecules reach thermal equilibrium approximately 4 ps after striking the surface at tphase ¼ 1 ps. Furthermore, the absence of a non-thermal population in the distribution of molecular speed indicates that all of the incident water molecules equilibrate to the liquid temperature. Therefore, S for water vapor on liquid water is determined to be 1.0 at 300 K. This value is in agreement with the value reported from the droplet train apparatus experiment  and other experimental work . The partitioning of the total kinetic energy during the thermal accommodation process is determined by resolving Eq. (4) into a sum of functions for each component of the kinetic energy
d KEtot ðtÞ ¼ dKExtr ðtÞ þ dKEytr ðtÞ þ dKEztr ðtÞ þ dKErot ðtÞ:
In Eq. (5), ‘xtr’, ‘ytr’ and ‘ztr’ represent the x, y and z translational kinetic energy, respectively, and ‘rot’ represents the rotational kinetic energy. Each term in Eq. (5) is similar in form to Eq. (4), except that the equilibrium ensemble average for the rotational kinetic energy and each component of the translational kinetic energy are 32kB T and 12kB T , respectively. The contributions to dKEtot ðtÞ from the components in Eq. (5) are shown in Fig. 4 for the adsorbed trajectories. The initial rise in the total kinetic energy is attributed to the translational kinetic energy in the direction normal to the liquid surface (z). The molecule is accelerated towards the surface as it approaches. The translational kinetic energies in the directions parallel to the surface (x and y) are not excited. There is a minor excitation of the rotational kinetic energy that is correlated with the incoming molecule striking the surface at tphase ¼ 1 ps. 3.3. Mass accommodation From the 250 initial conditions, two of the trajectories directly scatter. While one may immediately assume that the value of a is 0.99 based on this result, the long time behavior of the adsorbed trajectories that have not entered the bulk liquid within 90 ps must be considered since desorption into the gas phase may occur on a longer time scale than simulated. The correction to a based on desorption into the gas phase is the same as that employed by Wilson and Pohorille  for ethanol on water. This correction, p, is given by
Fig. 4. The resolved components of dKEtot ðtÞ as a function of time for the average of the 248 trajectories that resulted in adsorption. In each panel, the component (solid line) shown is labeled at the top and the total (dotted line) is shown for comparison. The time axis refers to the time passed since the incident water molecule reached below 21 A.
ebDAsurf!air ; ð6Þ ebDAsurf!bulk 1 where b ¼ ðkB T Þ , DAsurf!air is the free energy diﬀerence for transfer from the liquid surface to air and DAsurf!bulk is the free energy diﬀerence for transfer from the liquid surface to the bulk liquid. These free energy diﬀerences are obtained from the free energy proﬁle for the transfer of a water molecule across the air/liquid interface of water. While the free energy proﬁle for the POL3 water model is not available to our knowledge, the free energy proﬁles for a non-polarizable  and a polarizable  water model are available. Both proﬁles exhibit good agreement with the experimental solvation free energy, DAair!bulk , value of )6.3 kcal/mol at 298 K  and predict that a solute water molecule does not preferentially bind at the surface of liquid water. Since the free energy proﬁle monotonically decreases as the water molecule approaches the bulk liquid from air, DAsurf!bulk is taken to be equal to zero and the denominator in Eq. (6) is unity. If DAsurf!air is taken to be the experimental value of )6.3 kcal/mol, p is on the
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255
order of 105 . This minimal correction to a indicates that desorption of the adsorbed molecules into the gas phase from the surface is negligible and that a is indeed 0.99 for water vapor on liquid water at 300 K. This result is in agreement with a recent MD study of water mass accommodation by Morita et al.  using the SPC/E water model, which reports a > 0:99 at 273 K. Also using the SPC/E water model, an MD study by Tsuruta and Nagayama  determines a ¼ 0:961 at 330 K. The ﬁndings of those simulations and the results presented in this Letter indicate that a for water vapor on liquid water is temperature independent below 300 K and the negative temperature dependence begins in the range from 300 to 330 K. This is in clear disagreement with the temperature dependence of a determined by the droplet train apparatus experiment . In order to address this discrepancy, the MD simulation result by Morita et al.  is coupled with a ﬂuid dynamics simulation of the droplet train ﬂow tube. The results from the ﬂuid dynamics simulation indicate that the experimentally determined a is in the range from 0.2 to 1.0 at 273 K. The lack of agreement between the droplet train apparatus experiments and MD simulations is not exclusive to water vapor on liquid water. Determinations of a for ethanol [10,14] and methanol  on liquid water from MD simulations are near unity, while droplet train apparatus experiments predict a to be 0.056 and 0.049 for methanol and ethanol, respectively, at 273 K . The discrepancy between the droplet train experiments and MD simulations may be attributed to the lack of a correct description in the simulations for the intermolecular interactions of the species being considered. However, in the case of water, both models that are used, POL3  and SPC/E , exhibit good quantitative agreement with experimentally determined thermodynamic and structural properties of water [17,28,29]. Furthermore, diﬀerent simulation results for a are consistent, regardless of which water model is used [8,9]. It is unlikely that the empirical description of intermolecular interactions, which is a standard procedure in classical MD simulations, is the source of disagreement between determinations of a from experiments and simulations. The ﬁnal aspect of the mass accommodation process considered is the time scale for the incident molecule to obtain the equilibrium solvent structure after adsorbing to the liquid surface. To obtain information about this time scale, the oxygen–oxygen radial distribution function is calculated for the bulk and interface water molecules from the 250 ps equilibrium simulation (data not shown). The radial distribution functions are integrated to obtain the equilibrium number of water molecules in the ﬁrst, hn1 i, and second, hn2 i, solvent shells for a water molecule located in bulk water and at the air/water interface. These values are shown in Table 1.
Table 1 Equilibrium number of water molecules in the ﬁrst and second solvent shells of a water molecule in the bulk and at the air/water interface
Fig. 5 shows the time-dependent number of water molecules in the ﬁrst, n1 ðtÞ, and second, n2 ðtÞ, solvent shells for the 248 incident water molecules that adsorbed. Upon striking the surface at tphase ¼ 1:0 ps, the incident water molecules have 1.4 water molecules in the ﬁrst solvent shell, indicating that there are either one or two water molecules participating in the initial impact event, in addition to the incident molecule. The equilibrium values from Table 1 for a water molecule located at the air/water interface are obtained approximately 4 ps after striking the surface. The asymptotic values for both solvent shells are intermediate to the bulk and interface values shown in Table 1 since the average position of the incident water molecules corresponds to a water density that is approximately 90% of the bulk value (Fig. 3).
tphase (ps) Fig. 5. The time-dependent number of water molecules in the ﬁrst (thick line) and second (thin line) solvent shells of the incident water molecules averaged over the 248 trajectories that resulted in adsorption. The time axis refers to the time passed since the incident water molecules reached below 21 A.
J. Vieceli et al. / Chemical Physics Letters 393 (2004) 249–255
4. Conclusions Molecular dynamics computer simulations are used to determine the thermal and mass accommodation coeﬃcients for water vapor on liquid water at 300 K. S and a are determined to be 1.0 and 0.99, respectively. In addition, signiﬁcant molecular insight into the uptake process of a gas phase molecule with thermal impact velocity at a liquid surface is gained on a picosecond time scale. The process is characterized by an increase in the translational kinetic energy in the direction normal to the liquid surface that reaches a maximum upon impact with the surface. This collision induces a slight rotational excitation and the total excess kinetic energy is dissipated within 4 ps of striking the surface. Furthermore, the equilibrium number of water molecules in the ﬁrst and second solvent shells is obtained 4 ps after striking the surface. In an eﬀort to make a general conclusion regarding mass accommodation coeﬃcients, MD simulations have predicted that a is near unity for several diﬀerent gas phase species impinging on liquid water at a temperature near 300 K. The gas phase species that have been studied are water [8,9], ethanol [10,14], methanol , HO2  and ethylene glycol . This suggests that the mass accommodation coeﬃcient is not sensitive to speciﬁc features of the interaction between the incident molecule and water, such as the number of hydrogen bonds that the incident molecule can form or the solvation energy. A determination of a ¼ 0:26 for the OH radical on liquid water at 300 K from an MD simulation is in disagreement with the general result above . However, the OH radical force ﬁeld in that study underestimates the solvation free energy and the OH–H2 O dimer energy. Preliminary results with a revised OH force ﬁeld indicate that a is indeed closer to unity  and highlight the importance of accurate molecular force ﬁelds for MD simulations of the accommodation process.
Acknowledgements We thank Dr. Douglas R. Worsnop for encouraging us to perform this study. A Collaborative Research in
Chemistry grant from the National Science Foundation (CHE-0209719) funds this work. References  P.Y. Chuang, R.J. Charlson, J.H. Seinfeld, Nature 390 (1997) 594.  A. Nenes, S. Ghan, H. Abdul-Razzak, P.Y. Chuang, J.H. Seinfeld, Tellus 53B (2001) 133.  C.E. Kolb, P. Davidovits, J.T. Jayne, Q. Shi, D.R. Worsnop, Prog. React. Kinet. Mech. 27 (2002) 1.  R.A. Shaw, D. Lamb, J. Chem. Phys. 111 (1999) 10659.  G. Sageev, R.C. Flagan, J.H. Seinfeld, S. Arnold, J. Colloid Interface Sci. 113 (1986) 421.  Y.Q. Li, P. Davidovits, Q. Shi, J.T. Jayne, C.E. Kolb, D.R. Worsnop, J. Phys. Chem. A 105 (2001) 10627.  M. Mozurkewich, Aerosol Sci. Technol. 5 (1986) 223.  A. Morita, M. Sugiyama, H. Kameda, S. Koda, D.R. Hanson, J. Phys. Chem. B (2004) in press.  T. Tsuruta, G. Nagayama, J. Phys. Chem. B 108 (2004) 1736.  M.A. Wilson, A. Pohorille, J. Phys. Chem. B 101 (1997) 3130.  A. Morita, Chem. Phys. Lett. 375 (2003) 1.  A. Morita, Y. Kanaya, J.S. Francisco, J. Geophys. Res. 109 (2004) D09201.  M. Roeselova, P. Jungwirth, D.J. Tobias, R.B. Gerber, J. Phys. Chem. B 107 (2003) 12690.  R.S. Taylor, B.C. Garrett, J. Phys. Chem. B 103 (1999) 844.  L.X. Dang, B.C. Garrett, Chem. Phys. Lett. 385 (2004) 309.  P. Jungwirth, D.J. Tobias, J. Phys. Chem. B 106 (2002) 6361.  J.W. Caldwell, P.A. Kollman, J. Phys. Chem. 99 (1995) 6208.  D.R. Lide, CRC Handbook of Chemistry and Physics, CRC Press, Boca Raton, FL, 1992.  L.X. Dang, T. Chang, J. Chem. Phys. 106 (1997) 8149.  Q. Du, R. Superﬁne, E. Freysz, Y.R. Shen, Phys. Rev. Lett. 70 (1993) 2313.  F.J. Vesely, J. Comput. Phys. 24 (1977) 361.  P. Ahlstrom, A. Wallqvist, S. Engstrom, B. Jonsson, Mol. Phys. 68 (1989) 563.  J.P. Ryckaert, G. Ciccotti, H.C. Berendsen, J. Comput. Chem. 23 (1977) 327.  U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577.  D.A. Pearlman et al., Comput. Phys. Commun. 91 (1995) 1.  A. Ben-Naim, Y. Marcus, J. Chem. Phys. 81 (1984) 2016.  H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269.  J. Alejandre, D.J. Tildesley, G.A. Chapela, J. Chem. Phys. 102 (1995) 4574.  R.S. Taylor, L.X. Dang, B.C. Garrett, J. Phys. Chem. 100 (1996) 11720.  M. Roeselova, J. Vieceli, L.X. Dang, D.J. Tobias, (2004) in preparation.