An Internally Consistent Method for the Molecular Dynamics

482

J. Phys. Chem. B 2009, 113, 482–486

An Internally Consistent Method for the Molecular Dynamics Simulation of the Surface Tension: Application to Some TIP4P-Type Models of Water Raymond D. Mountain† Chemical and Biochemical Reference Data DiVision, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8320 ReceiVed: February 12, 2008; ReVised Manuscript ReceiVed: October 28, 2008

A set of molecular dynamics simulations of three members of the TIP4P model family for water (TIP4P, TIP4P/2005, TIP4P-i) are performed to illustrate the appropriate procedures to employ when simulating the planar interface between coexisting liquid and vapor phases. A method for efficiently including the longrange parts of the intermolecular dispersion interactions introduced by Janecˆek (J. Phys. Chem. B 2006, 110, 6264) is used. This method produces an interfacial region that is internally consistent with the dynamics of the molecules. A set of procedures are presented to verify that the surface tension estimate has converged and to obtain the underlying uncertainty of the surface tension estimate. The surface tension for the models is determined for temperatures between 300 and 550 K with an uncertainty of about 2 mN/m for this temperature range. The uncertainty is not sensitive to the temperature, but the relative uncertainty increases strongly with temperature. The surface tension is not sensitive to the system size in this temperature range provided at least 500 molecules are used, and it is shown that a thermostat does not bias the computed values. I. Introduction This paper contains a discussion of how molecular dynamics simulations of molecular, liquid-vapor planar interfacial systems can be performed in an internally consistent, computationally efficient manner with accurate inclusion of long-range interactions. This work extends the discussion of problems associated with the simulation of liquid-vapor interfaces recently published.1,2 The systems used to illustrate the method are from the TIP4P family of model potentials, namely the TIP4P,3 the TIP4P/2005,4 and the TIP4P-i5 water models, with emphasis on determining the surface tension and the density profile. The TIP4P/2005 model has been reported to closely reproduce the experimental values of the surface tension of water between 300 and 550 K6,7 so it is appropriate to use it and closely related models to illustrate the objective of this paper. The main purpose is to clarify how to make interfacial simulations in a way that the resulting density profile is consistent with the long-range dispersion and Couloumb interactions used in the simulation. Molecular dynamics simulations of the liquid-vapor interface have been used extensively to determine the surface tension of the interface.1 Most simulations have used a truncated potential in the simulation that determines the density profile and the surface tension. In many cases, some sort of long-range interaction correction was made to the surface tension after the simulation was completed with the intent of determining the surface tension for the full (rather than the truncated) model potential.2 This approach suffers from a serious flaw: the density profile obtained in the simulation is not consistent with the profile for the full model potential since it was determined with the truncated one. The correction applied at the end of the simulation uses the density profile for the truncated interaction and therefore has the same inconsistency problem. The inconsistency occurs because the forces acting on a particle from †

E-mail: [email protected].

10.1021/jp8012514

distant neighbors do not cancel out if the fluid is spatially inhomogeneous. These forces do cancel for a homogeneous fluid. Recently, Janecˆek developed a method for removing this inconsistency between the density profile and the surface tension.8 It involves making an approximate estimate of the contribution of the long-range part of the dispersion interaction to the dynamics of the system during the simulation. The direct interactions are truncated at some distance RC that is several times the effective size of the interacting particles. The simulation cell is divided into slabs oriented parallel to the liquid-vapor interface with thickness small compared to the size of the particles. During each step of the simulation, the density of particles in each of the slabs is determined and used, as described in the next section, to determine the forces on each of the atoms due to the long-range dispersion interactions. The density information is also used to determine the long-range contributions to the energy and the pressure tensor components. This results of this method have been extensively compared with surface tension estimates for the Lennard-Jones fluid obtained by other methods. The agreement is good, well within the uncertainties of the methods.9 The tests were applied to the one-component Lennard-Jones model, but the extension to mixtures is straightforward and is not discussed here. The development of approximate corrections for mixtures is a complex problem10 and is completely avoided using the Janecˆek approach. In this article, we apply this method to water models in the TIP4P family where the nonbonded, intermolecular interactions have the Lennard-Jones plus Coulomb form. The long-range part of the Coulomb interactions is treated using the Ewald summation method as this is a well-established procedure. The use of the Janecˆek method to include the long-range part of the Lennard-Jones interaction is new and thereby eliminates the internal inconsistency between the calculated density profile, the surface tension, and the model used to generate them.

This article not subject to U.S. Copyright. Published 2009 by the American Chemical Society Published on Web 12/16/2008

TIP4P-Type Models of Water The Janecˆek method is not the only possible way to account for the long-range part of the interactions. There is an Ewald summation-based method that in principle accounts for the longrange parts.11 This method is computationally much more expensive than the Janecˆek method and yields equivalent results for the surface tension if the appropriate number of wave vectors is used.9 Convergence of the summation method requires significantly more wave vectors than are needed for the Coulomb interactions so the computational effort is large and not commensurate with the effort needed. There have been a number of simulation studies of the surface tension of water for several models showing that the commonly used water models underestimate the surface tension by 15-30% of the experimental value.1,6 The TIP4P/2005 model is an exception. The liquid and vapor coexisting densities for the model are reported to be in good agreement with experimental values12 although there is a conflicting report in the literature5 that will be discussed below. Procedures for determining the adequacy of the sampling of configurations are considered along with the corresponding estimates for the uncertainty in the computed value of the surface tension. The influence of both system size and thermostats on the surface tension estimates is examined briefly. The methods described here are also applicable to Monte Carlo simulations of planar interfacial systems. The simulation methods and a discussion of the uncertainty in the simulation values of the surface tension are contained in section II. The states considered here are described in section III, and the values of the surface tension for the NVT and NVE ensembles are compared. The results of this study are discussed in section IV and are compared with prior results. Finally, some conclusions are presented about how molecular dynamics and Monte Carlo simulations that generate the interface explictly should be performed. II. Simulation Methods and Uncertainty Estimates The interaction between water molecules is described by the planar TIP4P-type models. These are rigid, four-site models with the OH distance of 0.09572 nm, and the HOH angle is 104.52°. There is also a fourth site, the “M” site, located on the bisector of the HOH angle between the oxygen and hydrogen sites. There is a Lennard-Jones interaction between the oxygen sites and Coulomb interactions between the hydrogen and M sites. There are four parameters that define the various TIP4P-type models. These are the Lennard-Jones
and σ parameters, the charge, qH, on the hydrogen sites, and the distance, δM, between the oxygen site and the M site. The change on the M site is -2 times the charge on a hydrogen site. For the TIP4P model,3 the parameters are
) 0.649 kJ/mol, σ ) 0.3154 nm, qH ) 0.52e, and δM ) 0.015 nm. The TIP4P/ 2005 model parameters are 0.7743 kJ/mol, 0.31589 nm, 0.5564e, and 0.01546 nm, respectively.4 The TIP4P-i parameters are 0.7801 kJ/mol, 0.3163 nm, 0.5673e, and 0.01556 nm.5 Simulations were made for each of the models with the emphasis on the TIP4P/2005 model as this model is reported6 to produce surface tension values in good agreement with the experimental values. The results for the TIP4P model are included for completeness. Some simulations were also performed using the TIP4P-i model.5 It was developed to provide a more accurate representation of the liquid-vapor equilibrium envelope. These authors reported a critical temperature for the TIP4P/2005 model of 617 K5 rather than 640 K as reported by the developers of that model.12 The critical temperature of the TIP4P-i model is 639

J. Phys. Chem. B, Vol. 113, No. 2, 2009 483 K.5 A comparison of the temperature variation of the surface tension for these two models will be made in sections III and IV. The specific states examined are listed in section III. Most of the states simulated consist of 500 water molecules in a rectangular simulation cell with dimensions LX ) LY and LZ ) 3LX. Periodic boundary conditions are applied in all directions. Several cases with 1000 molecules and LX and LZ appropriately increased to maintain the same global density were examined. With a few exceptions, the temperatures of the systems were maintained by separate Nose´-Hoover thermostats for the translational and rotational degrees of freedom13 (NVT ensemble). The time constants for the thermostats were set to ∼0.5 ps as this was found to provide good temperature control. A few constant energy simulations (NVE ensemble) were made to examine the sensitivity of the simulation value of the surface tension to the use of thermostats. The equations of motion of the molecules were integrated using a velocity-Verlet type algorithm that uses quaternions14-17 to describe the orientation of the molecules.18 For constant temperature simulations, the time step used was 1 fs, and for constant energy simulations, the time step was 0.5 fs. The shorter time step for the NVE simulations is needed to maintain energy conservation at the elevated temperatures simulated here. The direct interactions are truncated at RC ) 0.5LX. The longrange parts of the Coulomb interactions are evaluated using the Ewald method as discussed by Alejandre et al. and implemented by Wemhoff.19,20 The long-range part of the Lennard-Jones interactions is evaluated using the method of Janecek8 with slabs of thickness δz ) 0.01 nm. The contributions to the energy, force, and pressure tensor components associated with a molecule due to the long-range dispersion interactions are determined by taking the volume integral of the energy, etc., over each slab. The distance in the z-direction from the molecule to the slab is held fixed, and only that portion of the slab that is more distant than RC from the molecule is included in the integration. The volume integrals are weighted by the number density of molecules in the slab to obtain the energy, etc., for the molecule in question. The total effort involves sums over all molecules and all slabs. The actual volume intergrals are known analytically for the Lennard-Jones potential and are derived in the paper where the method is developed8 and are not repeated here. Note that periodic boundary conditions in the z-direction are used when determining the z-distance from the molecule to the slab. There are two approximations involved in this process. The first is the assumption that the molecules in the slab are uniformly distributed, and the second is that the contributions from slabs further from a molecule than 1/2 the long dimension of the simulation cell can be neglected. Given that this method yields results in good agreement with other methods indicates that these approximations do not cause problems. The long-range Lennard-Jones contributions to the surface tension are important for the system sizes simulated and must be included if reliable values for the surface tension are to be realized. In the method used here, the long-range terms are included in the forces, pressure tensor components, and energy at each time step of the simulation. It is not appropriate to fit a functional form, such as a hyperbolic tangent, to the density profile and use that fit profile to make a mean-field adjustment to the pressure tensor components after the simulation is completed.21 Doing so introduces an internal inconsistency between the intermolecular interactions and the resulting density profile. By including the long-range part of the interactions in

484 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Figure 1. Average mass density profile for TIP4P/2005 water at 450 K is shown as a function of z, the direction normal to the interface. The solid line is the simulation result, and the circles are a tanh fit to the solid line.

Mountain

Figure 2. Cumulative time average surface tension for a 500-molecule system of the TIP4P/2005 model at 300 K is shown for three successive runs of 1.5 ns duration each. There is no offset of the curves. The solid line is for the first run, the dashed line is for the second run, and the dash-dot line is for the third run in the sequence.

the forces and pressure tensor components, the density profile and the computed surface tension are internally consistent with the dynamics of the molecules. An estimate of the average liquid density, Fl, can be obtained by fitting the tanh profile form

F(z) ) 0.5(Fl + Fg) + 0.5(Fl - Fg) tanh((z - z0)/d) to the average density profile where Fg is the gas density, z0 is the midpoint of the interface, and d specifies the width of the interface. A good fit to the simulation profile using this form is not a good indicator of the gas density when it is small. The 500-molecule simulations are run for 1.5 ns once a stable two-phase system is established. As is shown below, runs of this duration are the minimum length run needed if stable, reasonably accurate estimates for the surface tension are to be realized for a 500-molecule system. The surface tension γ is determined by evaluating the Kirkwood-Buff result22

γ ) [〈pzz〉 - 0.5(〈pxx〉 + 〈pyy〉)]/(2A) where A ) LXLY is the cross-sectional area of the cell and 〈paa〉 is the average aa-component of the pressure tensor. The extra factor of 2 is due to the presence of two liquid-vapor interfaces. The pressure tensor components are evaluated using the molecular virial.23 The fluctuations in the pressure tensor components are large, and long runs are needed to adequately sample the fluctuations. Stable results were obtained for a 1000molecule system with runs of 1 ns duration. An example of the configuration of the molecules is shown in Figure 1 for the TIP4P/2005 model where the average density profile is shown for 450 K. The average is taken over slices normal to the z-axis of 0.1 nm thickness and over the duration of the simulation. The circles represent a tanh fit to the density profile. (The right-hand part of the fit subtracts from rather than adds to the tanh term with the appropriate value for z0.) During a run, the pressure tensor components are averaged in blocks of 1 ps duration for later processing. When a run is completed, the surface tension, averaged over the individual block values, is determined as well as the standard error24 of the average of the block values. The standard error for a 500molecule system is ∼1.5 mN/m for a 1.5 ns run over the temperature range investigated. This provides one estimate of the uncertainty in the calculated value of the surface tension. A second estimate of the uncertainty is obtained by performing three sequential runs. The spread in the average surface tension values is on the order of 2-3 mN/m or less. A third estimate of the uncertainty can be obtained by examining the cumulative time average of the 1 ps values as the system evolves in time.

Figure 3. Distribution of 1 ps block averages of the surface tension is shown for the three successive runs in Figure 2. The sequence is the same as in Figure 2; the lowest one corresponds to the lowest curve in Figure 2. The distributions are offset by 0.01 and 0.02 from the first one.

Figure 4. Time correlation function of 1 ps block averages of the surface tension is shown for one of the 350 K TIP4P/2005 cases.

This is illustrated for three runs at 300 K in Figure 2. Good convergence of the estimate of the surface tension is indicated by the length of the approximately flat portion of γ(t), the cumulative average of γ at time t, for longer times. To gain some insight into why these three cases appear to converge to slightly different values, we examine the distribution of block averages that are used to construct γ(t). The distribution of the 1 ps block averages, P(γ), for the three time averages shown in Figure 2 is displayed in Figure 3. The distributions are offset by 0.01 and 0.02, respectively. The shape of P(γ) illustrates the earlier comment that the fluctuations in the pressure tensor components are large. The distributions for the cases in Figure 2 were generated in different sequences with

TIP4P-Type Models of Water

J. Phys. Chem. B, Vol. 113, No. 2, 2009 485

TABLE 1: 500-Molecule States Simulated in the NVT Ensemble with the Cell Dimension and the Surface Tension Estimate Listeda T, K

LX, nm

γ, mN/m

γ, mN/m

300 300 300 300 350 350 400 400 450 450 450 450 500 500 550 550 600

2.52 2.52 2.52 2.643 2.52 2.643 2.52 2.643 2.52 2.643 2.643 2.643 2.52 2.643 2.52 2.643 3.175

69.3 (2.8) 71.2 (2.8) 68.3 (3.0) 67.7 (2.4) 63.0 (2.8) 61.7 (2.4) 51.8 (2.8) 52.3 (2.6) 42.1 (2.8) 41.9 (2.4) 42.3 (2.4) 42.4 (2.6) 31.7 (2.8) 31.7 (2.4) 19.5 (2.4) 20.4 (2.0) 9.1 (2.2)

74.1 (2.8)

γ, nm/m

T, K

LX, nm

γ, mN/m

57.1 (2.2)

301 403 497 499

2.643 2.643 2.52 2.52

68.0 (3.2) 50.1 (3.0) 33.4 (4.6) 31.1 (3.2)

65.4 (2.6) 47.8 (2.2) 57.1 (3.0) 37.4 (2.2)

a Two times the standard error is shown in parentheses. These simulations were made to examine the sensitivity of the results to the thermostats used in the simulations listed in Table 1.

46.7 (2.8) 27.3 (2.2) 36.0 (2.8) 15.8 (1.8) 25.6 (2.6) 13.6 (2.6)

a Two times the standard error is shown in parentheses. The states run more than once are listed separately to illustrate the variability of surface tension estimates for 1.5 ns duration simulations. Column 3 contains the surface tensions for the TIP4P/ 2005 model, column 4 contains the surface tensions for the TIP4P-i model, and column 5 contains the surface tensions for the TIP4P model. The 600 K states are for 1000-molecule systems of 1.0 ns duration.

the lowest curve case predominantely sampling the larger values earlier in the run and then the smaller values later in the run while the other two cases sampled the distribution more uniformly over the course of the run. In all three cases, the distribution of 1 ps block averages is not a smooth curve, indicating that the sampling of the fluctuations in the pressure components is not complete. The standard error estimate assumes that the 1 ps duration block averages are statistically independent. This can be checked by constructing the time correlation function of the block averages, C(τ)25 N

C(τ) )

TABLE 2: 500-Molecule States Simulated in the NVE Ensemble with the Cell Dimension and the Surface Tension Estimate Listed for the TIP4P/2005 Modela

∑

1 (γ(ti) - 〈γ〉)(γ(ti + τ) - 〈γ〉) N i)1

where the sum is over the time origins, γ(ti) is the block average for time ti, and 〈γ〉 is the average over the blocks. When C(τ) approaches zero and remains small, the blocks with that τ and larger separation are statistically independent. An example of C(τ) is shown in Figure 4 for 350 K. For the various cases examined here, C(τ) becomes small when τ is 2 to 3 ps or larger, except for 300 K where τ is about 5 ps. This indicates that the standard error underestimates the uncertainty in the computed values of the surface tension by about a factor of 2. Given the spread in final values of the surface tension in Figure 2, it is evident that the standard error underestimates the uncertainty in the simulation result, at least for runs of 1.5 ns duration. A better uncertainty estimate would be twice the standard error as that would be consistent with the spread of averages in Figure 2 and with the results of the time correlation function analysis. III. Results The states of TIP4P/2005 water examined are listed in Table 1 (NVT ensemble) and Table 2 (NVE ensemble) for 500molecule systems. The simulation results for the TIP4P and the

Figure 5. Surface tension results for the TIP4P/2005 model and for the TIP4P-i model are compared with the experimental values (solid line). The dashed line and squares are for TIP4P/2005 with LX ) 2.52 nm, and the dash-dot line is for that model with LX ) 2.643 nm. The triangles and the dash-dash-dot line are for the TIP4P-i model.

TIP4P-i models are also listed in Table 1. In both tables the uncertanties listed in parentheses are 2 times the standard error based on block averages as discussed in section II. The states are specified in terms of the temperature and the X-dimension of the simulation cell. From Table 1 we see that a 10% change in the interfacial area does not result in variations in the surface tension larger than the underlying uncertainty in the individual estimates. This is evident in Figure 5 where the simulation results are compared with the experimental values of the surface tension of water.26 The TIP4P/2005 surface tension values are quite close to the experimental values while the TIP4P-i values are systematically larger than the experimental ones and the TIP4P values (not shown) are smaller than the experimental ones. In addition, several simulations of a 1000-molecule system with LX ) 3.175 nm were run for 1 ns. The 1000-molecule results are consistent with the 500-molecule results so 500 appears to be an adequate size for temperatures up to 550 K. This is not the case for higher temperatures. A 500-molecule case was run at 600 K for the TIP4P/2005 model, and the liquid phase turned out to be poorly defined, unlike the density profile shown in Figure 1. Instead of uniform density regions for the liquid and vapor phases, there was a broad maximum so that the “interface” is ill-defined. The critical temperature for the TIP4P/2005 model is 640 K.12 For systems within about 10% or less of the critical point, a larger number of molecules or much longer simulation times (or both) are needed to obtain reliable results. A 1000-molecule system was run at 600 K, and reasonably well-defined liquid and vapor regions were realized. The 600 K surface tensions are listed in Table 1. By comparing the results in Table 2 with the results for the equivalent states in Table 1, we see that the use of a Nose´Hoover thermostat does not shift the calculated values of the surface tension from the values obtained without a thermostat.

486 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Figure 6. Surface tensions are shown as functions of (1 - T/Tc)1.256 where the critical temperatures used are discussed in the text. The solid line is the experimental curve, The dashed line with circles is the TIP4P curve, the dash-dot line with squares is the TIP4P/2005 curve, and the dash-dash-dot line with triangles is the TIP4P-i curve. The error bars for the simulation results are taken from Table 1.

IV. Discussion The published values of the surface tension of TIP4P/2005 water6,7 are consistent with the values obtained here. Within the estimated uncertainties of 2-3 mN/m in the simulation results, there is good agreement between the simulation and experimental values as shown in Figure 5. There is an alternative way to examine the temperature dependence of the surface tension. At elevated temperatures, the surface tension of water varies as (1 - T/Tc)1.256, where Tc is the critical temperature.26 This representation of the surface tension is displayed in Figure 6 for the experimental values, for the TIP4P values (Tc ) 588 K), for the TIP4P/2005 values with 640 K for the critical temperature, and for the TIP4P-i model. The experimental and the TIP4P/2005 with Tc ) 640 K12 surface tensions extrapolate smoothly to zero at the critical temperature, as they must, while the TIP4P-i surface tension with Tc ) 639 K appears to extrapolate to a finite value at Tc. The 640 K value of the critical temperature of the TIP4P/2005 model is consistent with the temperature variation expected for the surface tension while with the lower critical temperature mentioned previously5 is not. This suggests that the reported critical temperature of a model can be tested for consistency with the surface tension and that Tc for the TIP4P-i model is too low. This study indicates that the following conditions should be satisfied if reliable liquid-vapor interface simulations are made. The system should contain at least 500 molecules, and the width of the simulation cell should be on the order of 8 times the Lennard-Jones sigma. The length should be 3 times the width. It is necessary to perform simulations of at least 1.5 ns duration, and the convergence of the surface tension estimate should be checked using analyses of the sort that produced Figures 2 and

Mountain 3. It is likely that significantly longer runs will be required for larger molecules.27 The uncertainty in the estimate of the surface tension should be determined using block average values to determine the standard error, and C(τ), the time correlation function of the block averages, should be used to scale the standard error as appropriate. It is permissible to use a Nose´-Hoover thermostat to control the temperature in a molecular dynamics simulation. It is necessary to accurately include the long-range parts of the intermolecular interactions. The Janecˆek method is the appropriate one for the Lennard-Jones interaction, and for the Coulomb interactions the Ewald method (or its equivalent) with an adequate number of wave vectors is the appropriate method. These treatments of long-range interactions produce results that are internally consistent with the underlying dynamics of the molecules. The high-temperature variation of γ as a function of (1 - T/Tc)1.256 can provide a consistency check for γ and Tc. References and Notes (1) Chen, F.; Smith, P. E. J. Chem. Phys. 2007, 126, 221101. (2) Duque, D.; Vega, L. F. J. Chem. Phys. 2004, 121, 8611. (3) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (4) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (5) Chialvo, A. A.; Barto´k, A.; Baranyi, A. J. Mol. Liq. 2006, 129, 120. (6) Vega, C.; de Miguel, E. J. Chem. Phys. 2007, 126, 154707. (7) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154716. (8) Janacˆek, J. J. Phys. Chem. B 2006, 110, 6264. (9) Shen, V. K.; Mountain, R. D.; Errington, J. R. J. Phys. Chem. B 2007, 111, 6198. (10) Chen, F.; Smith, P. E. J. Phys. Chem. B 2008, 112, 8975. (11) Karasawa, N.; Goddard, W. A. J. Phys. Chem. 1989, 93, 7320. (12) Vega, C.; Abascal, J. L. F.; Nezbeda, I. J. Chem. Phys. 2006, 125, 034503. (13) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (14) Evans, D. J.; Murad, S. Mol. Phys. 1977, 34, 327. (15) Mackay, A. L. Acta Crystallogr. 1984, A40, 165. (16) Sonnenschein, R. J. Comput. Phys. 1985, 59, 347. (17) Rapaport, D. C. J. Comput. Phys. 1985, 60, 306. (18) Martys, N. S.; Mountain, R. D. Phys. ReV. E 1999, 59, 3733. (19) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 2002, 102, 4574. (20) Wemhoff, A. P.; Carey, V. P. Int. J. Thermophys. 2006, 27, 413. (21) in ’t Veld, P. J.; Ismail, A. E.; Grest, G. S. J. Chem. Phys. 2007, 127, 144711. (22) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1949, 17, 338. (23) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987; pp 46-49. (24) Hogg, R. V.; McKean, J. W.; Craig, A. T. Introduction to Mathematical Statistics, 6th ed.; Pearson Prentice Hall: Upper Saddle River, NJ, 2005. (25) Huang, D. M.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704. (26) White, H. J., Jr., Sengers, J. V., Neumann, D. B., Bellows, J. C., Eds.; Physical Chemistry of Aqueous Systems: Meeting the Needs of Industry; Proceedings of the 12th International Conference on the Properties of Water and Steam; Begell House: Redding, CT, 1995; pp A139-A142. (27) Ismail, A. E.; Tsige, M.; in ’T. Veld, P. J.; Grest, G. S. Mol. Phys. 2007, 105, 3155.

JP8012514