Faraday Symp. Chem. Soc. 1982 17,55-67 Thermodynamics of Cavity Formation in Water A Molecular Dynamics Study BY JOHAN P. M. POSTMA J. C. BERENDSEN HERMAN AND JAN R. HAAK Laboratory of Physical Chemistry The University of Groningen Nijenborgh 16 9747 AG Groningen The Netherlands Received 26th August 1982 Thermodynamic quantities related to the solvation of hydrophobic solutes in water can be approxi- mated by the application of scaled-particle theory. The crucial quantity is the Gibbs free energy of cavity formation. A series of six molecular-dynamics simulations of water including repulsive cavities of various sizes has been carried out. Using perturbation statistical mechanics the free energy has been derived as a function of cavity radius up to 0.32 nm (0.32 nm approach of water oxygens to the cavity centre).The free energy agrees well with predictions from scaled-particle theory and the experimental surface tension is predicted to within 5%. The radial distribution of water molecules with respect to the cavity has been determined for five cavity sizes; for one size (radius of approach of water oxygens of 0.3 nm) the orientational distribution and the residence-time distribution in the hydration shells has been determined. 1. INTRODUCTION The thermodynamics of hydrophobic hydration can be reasonably well understood by the application of scaled-particle theory (SPT). This theory originally devised for hard-sphere solutes and solvents,' has been applied to aqueous solutions by Pierotti 2*3 and refined for water by Stillinger.4 Several later applications and comments are a~ailable.~-~ SPT essentially computes the free energy of formation of a hard-sphere cavity of diameter a2 in a hard-sphere solvent of molecular diameter a and number density p by fitting a third-degree polynomial to the exact solution available for small cavity sizes.For application to real solutions a perturbation treatment is used to allow for realistic interactions between solute and solvent. The reasonable success of SPT for aqueous solutions of hydrophobic solutes suggests that the main factor determining the negative entropy of hydration is the packing of solvent molecules around the solute. The determining parameters are the hard-sphere size of the solvent set equal to the normal hydrogen-bonding approach of two water molecules and the experimental density of water.In section 2.1 the results of SPT are summarized. The applicability of SPT to water remains to be assessed because rather severe approximations are made. The most important of these is the hard-sphere assump- tion water may be expected to be capable of the formation of a variety of hydrogen-bonded structure (" icebergs '' in a rather dated terminology) that could well deviate in thermodynamic properties from hard-sphere solvent structures. In as far as the radial distribution function of liquid water gives a clue to such deviations corrections have been derived by Stilli~~ger.~ A second important approximation is that the introduction of solutes into hard-sphere cavities is treated as a perturbation of the hard- sphere cavity reference state.Since such a state does not contain water configurations protruding into the cavity soft potentials cannot be reliably introduced. Moreover SIMULATION OF CAVITY FORMATION solute perturbation is usually not small. In practice the entropic contribution to the perturbation is negle~ted.~ Molecular-dynamics (MD) simulations using a fast computer provide the means to study the process of cavity formation in detail. Both structural and dynamical properties of the hydration shell follow from the simulations. The treatment can incorporate any soft cavity potential. Most importantly if simulations are carried out over a range of solute sizes such that overlapping statistics result the full thermo- dynamic solvation quantities can be derived by integration over a scaling parameter starting at the pure liquid.Molecular dynamics is thus an extremely powerful way to answer any pertinent question on hydration and is limited in principle only by the accuracy of the potentials used. MD and Monte Carlo simulations on hydrophobic hydration have been reported.10-14 We have chosen the following approach. Using a simple but well tested point- charge model for water l5 we have carried out isothermal-isobaric simulations of 216 water molecules in a periodic box including repulsive cavities of five different sizes. The repulsive interaction has the form Vcw= L(B/r)'*,where L is a scaling parameter B is the same parameter as used in the repulsive term of the water-water potential and r is the distance from an oxygen to the cavity centre.Simulations extended over 30-35 ps. Details of the simulations are given in section 3. By choosing appropriate cavity sizes and by monitoring appropriate perturbation functions the free energy for the formation of a cavity of sizes up to 0.32 nm (radius of approach of water mole- cules) was obtained. The perturbation method by which these values were determined is described in section 2.3 while results are given in section 4.1. Structural properties of the cavity hydration are described in section 4.2 with emphasis on one simulation (with cavity approach radius of 0.3 nm). Of the various dynamic properties (trans- lational and rotational characteristics of the hydration layers) we analyse the residence- time distribution in the hydration shells of one cavity diameter in section 4.3.A full account of all structural and dynamical properties including a more accurate evalua- tion of thermodynamic properties over a more extended range of cavity sizes will be reported elsewhere. Finally in section 5 we discuss the results of these simulations in comparison with scaled-particle theory thereby assessing the limitations and accuracy of SPT as applied to water. 2. THEORETICAL 2.1. SCALED-PARTICLE THEORY lt3 The standard Gibbs free energy of transfer of one mole of a substance X from its ideal-gas phase at standard pressurep and temperature T to an infinitely dilute solution at the same temperature and pressure referred to a standard of unit mole fraction of X is given by3 AGO = AGc + AGi + RTln-RT POGl where AGc is the isobaric reversible work to create a cavity with given radius and at a fixed position of its centre AGi is the reversible work to charge the cavity with the real interaction between solute and solvent and G1 is the molar volume of the solvent.It is assumed that the internal partition function of the solute is the same in solution and in the gas phase. AGO relates to Henry's constant KH(the ratio between gas fugacity and mole fraction in a saturated dilute solution of the gas) AGO = RTln KH. J. P. M. POSTMA €3. J. C. BERENDSEN AND J. R. HAAK The reversible work to create a mole of hard-sphere cavities with radius of approach r is given by’ 4AG = -RT In (1 -y3p) (r < a/2) (3) and AG = KO+ K,r + K2r2+ K3r3 (r > a/2) (4) where 6 + 18(Y)I+npa2 1-Y 4 -np.(44 K3-3 Here a is the (hard-sphere) diameter of the solvent y = 7ca3p/6is the volume fraction of the solvent spheres p is the number density of the pure solvent and p the pressure. The terms involving p are negligible for cavities of molecular size at atmospheric pressure. The radius of approach r is the smallest distance between a solvent mole- cule and the cavity centre; it is related to the diameter a of a hard-sphere cavity by r = (a + a2)/2. Eqn (3) is exact (1 -4nr3p/3)is the probability that the cavity centre is not within a distance r of a solvent centre (r < a/2; r = a/2 corresponds to a point cavity with a2= 0) and this probability is equal to exp(-AGJRT) where AG is the reversible work to create such a cavity.The creation of a cavity also requires enthalpy because of the temperature dependence of p and y. Values of AG are given in fig. 1 for T = 305 K using a = 0.2875 nm and a density of 0.983 g CM-~(see section 4.1). 2.2. THERMODYNAMICS OF SOFT-CAVITY FORMATION In our simulations we have introduced a soft cavity with cavity-water interaction potential (see section 3) VCw= A(B/r)12 (5) where B = 0.3428 (kJ mol-1)1/12 nm. The scaling parameter ilcan also be regarded as a scaling factor for the size of the cavity. We define the thermalradius rth of the cavity as the radius at which Y, =kT rth= B(L/kT)1’12 (6) Just as in the SPT for hard spheres we can derive a rigorous expression for AG for small cavities.For this we use VCwas a perturbing potential energy term in the Hamil- tonian of pure water. As shown in subsection 2.3 G -Go = -RT In ( exp(-VcW/kT) (7) )o SIMULATION OF CAVITY FORMATION 't -10 t ti +* I 1 I 1 I 1 0.1 0.2 0.3 thermal radius/nm FIG.1 .-Free energy of cavity formation derived from simulations (A) compared with scaled-particle theory (B). Points C denote relative AG values obtained by perturbation of cavity diameter in each of the five simulations. Thermal radii of the simulations are indicated by arrows. where GIand Goare the free energies of cavity solution and pure solvent respectively and ( )o denotes an average over the isobaric-isothermal ensemble for pure water.Now for small cavities interaction is only appreciable at small distances r. As long as r < 0.5 rmin,where rminis the distance below which the radial distribution function of water shows no neighbours (rmin x 0.24 nm) interaction takes place with maximally one water molecule. If po(r) denotes the normalised probability distri- bution of the distance of the nearest water molecule to the (arbitrary) cavity centre we obtain W WP(-VCWlkT) ) 0 = I exP(- Vcw/kT)Po(r)dr. (8) 0 For small r,po(r)=4nr2p and for small R the exponential in the integrand changes from zero to one in this range. Taking into account that po(r) is normalized it follows that for small A W (exp(-Vc,/kT)) = 1 + p [exp(- Vc,/kT) -1]4nr2dr.(9) 0 The integral can be expressed in terms of a gamma function and is also known from computation of second virial coefficients for power-law repulsion.16 The result is J. P. M. POSTMA H. J. C. BERENDSEN AND J. R. HAAK where r(:) = 1.226 41. Except for a small numerical factor eqn (10) is equivalent to the hard-sphere expres- sion (3) if we rewrite eqn (10) with eqn (6) in terms of the thermal radius lim AGc = -RTln rth Expression (11) should be valid for thermal radii up to ca. 0.1 nm. For radii between 0.5 rminand (1/43) rmin(ca. 0.14 nm) extension is possible using the radial pair distribution curve g(r) of pure water;Ip4 above (1/43) rminmore than two water molecules can exist at the cavity radius and higher-order distributions are required.We derive values for AG by overlapping perturbations of MD simulations. 2.3. THERMODYNAMIC VALUES FROM PERTURBATIONS IN STATISTICAL REFERENCE ENSEMBLES Let us assume that we have a representative classical statistical ensemble (from MD or Monte Carlo simulations) available containing a sufficient number of states to yield accurate averages. We restrict ourselves to an isobaric-isothermal ensemble but the results are easily modified for other ensembles. Let the Hamiltonian be given by Xo(p 4). The free energy G is related to the isobaric partition function A G= -kTlnA (12) where Of course this integration cannot be carried out. Now consider a system with a perturbed hamiltonian Z1=Xo+ AV,where A V is a perturbation of the potential energy.The difference in free energy AG is AG = G1 -Go= -kT In (Al/Ao) (14) where A1 is similar to eqn (13) involving Z instead of Zo.It follows that or AG = -kT In ( exp(-V/kT) ) o. (16) If AV is the perturbation resulting from the introduction of or change in the proper- ties of one particle (cavity) AG applies to one particle; per mole of particles we have AG = -RTln ( exp(-AV/kT) ) o. (17) The averaging is over the unperturbed isobaric-isothermal reference ensemble. In practice it is possible to derive reliable values for AG if IAGl does not exceed 2kT; for larger deviations the reference ensemble contains an insufficient number of con- figurations that are representative for the perturbed ensemble. Values for the enthalpy H can be derived directly from each ensemble so that the SIMULATION OF CAVITY FORMATION need for perturbative treatment does not arise.One can easily derive however that The first term indicates the direct interaction and the second is a contribution from relaxation of the ensemble. 3. MODEL AND METHOD The model used for the water molecule is the SPC (simple point charge) model which is an effective pair potential derived from MD ~imulations.~~ It reproduces the radial distribution curve of liquid water with a second neighbour peak at the correct position but less pronounced than the X-ray results. Density is ca. 2% too low; energy is correct as well as specific heat and compressibility.18 The diffusion constant is ca. 50% too high.The model is similar to the TIPS mode1;19 it consists of a negative charge of -0.82e at the oxygen nucleus and two positive charges of 0.41e on tetrahedral positions at 0.1 nm from the oxygen. No dielectric constant is used. In addition the oxygens have a Lennard-Jones interaction Simulations were carried out on 216 molecules in a periodic cube with intermolecular interactions taken into account for all water molecules for which roo < 0.9 nm. Cartesian coordinates were used and the internal degrees of freedom were constrained by the SHAKE method.20*21 A time step of 0.001 ps was used which gave a (short time) fluctuation of the total energy better than 15% of the kinetic-energy fluctuation. An algorithm was used that implies a coupling of the system to a constant tem- perature and pressure bath.22 The algorithm is based on a leap-frog scheme in which at each step both velocity and coordinate scaling (and volume scaling) occurs.The scaling involves delayed coupling to constant temperature and pressure; the time constant of each coupling is adjustable. A long time constant means a weak coupling and slow restoration of temperature and pressure deviations but also a weak pertur- bation of the system. We used 0.4 ps for the coupling to constant temperature and 0.5ps for the coupling to constant pressure. Dynamic runs were initiated by gradually introducing a cavity starting from a pure water configuration in 1 or 2 ps and then equilibrating over 5 ps. Runs extended over 25 or 30 ps. Averages over the equili- brating 5 psdid not differ significantly from any 5 ps section of the run but were never- theless not used for obtaining final averages.The runs were subdivided into 5 or 6 sections of 5 ps each; averages Aiwere obtained for each subsection and final averages A' and their estimated errors AA were obtained from the N subaverages and their statistical spread IN A=-2 Ai (204 N i=l N The following potential between cavity centre and oxygen of water molecules J. P. M. POSTMA H. J. C. BERENDSEN AND J. R. HAAK was introduced for five values of A. The cut-off radius used for this interaction was also 0.9 nm. The cavity was kept stationary. 4. RESULTS 4.1. THERMODYNAMIC QUANTITIES Data for six MD runs including pure water are given in table 1.It is apparent that neither the volume nor the enthalpy shows significant dependence on cavity size. This is because the full fluctuations of the large system of 216 molecules are present in such overall quantities while the cavity concerns effects on a single-molecule scale. TABLERESULTS FROM MD RUNS A t-th length Ew-wPot Ec-wPot V H ~p0-5 Inm /ps /kJ mol-' /kJ mol-' /nm3 /kJ mol-' cm2s-I 0 0 25 -8983 (13) 0 6.57 (1) -7325 (15) 4.2 (6) 0.100 25 -8996 (10) 0.1 (1) 6.63 (2) -7345 (22) 5.9 (4) 10-3 0.178 30 -8997 (12) 1.3 (1) 6.62 (1) -7340 (13) 4.6 (6) 0.0325 0.238 30 -8992 (10) 2.9 (1) 6.61 (2) -7332 (11) 7.1 (12) 0.5 0.299 25 -8971 (13) 6.2 (9) 6.65 (5) -7306 (15) 6.8 (7) 1 0.317 25 -9013 (6) 6.9 (6) 6.65 (2) -7343 (9) 4.8 (3) a Temperature for all runs is 305 K (rt0.5 K).Pressures are between 30 and 50 bar not corrected for cut-off. Values in parentheses are r.m.s. errors in last digit determined from statistics of the averages over subsets of 5 ps each. Free energy of cavity formation was determined in the following way. For each value of rth the perturbation treatment of section 2.3 [eqn (17)] was carried out for several deviations of rth on both sides of the thermal radius. Thus overlapping extra- polations from adjacent simulations were obtained that enable construction of the AG curve over the full range of cavity sizes studied. For each perturbation the error in AG was determined according to eqn (20) by performing the averages separately over 5 ps sections. Fig. 1 shows the values of AG obtained in this way for the various simulations as well as the AG curve constructed from these.The error in AG accumulates when the number of successive overlap adjustments increase; it is indicated by error bars in fig. 1. The connection between the simulation at rth = 0.100 and pure water deserves extra attention since extrapolation from this simulation downward to water is not possible. This is due to the fact that a cavity does not contain any water molecules; hence statistics rapidly become poor when extrapolations to smaller cavity size are made configurations are then required that are not available in the reference simu- lation. This is true for all reference sizes and extrapolation upward is always more reliable then downward.Even the error gives an unreliable estimate when insufficient statistics exist. Thus for the connection between rth = 0.1 nm and water several more simulations would be required. Fortunately the theoretical value for AG of eqn (1 1) is reliable up to 0.1 nm and can be used to obtain the connection. This was in fact done; fig. 2 illustrates this point. Fig. 1 also gives values of AGc for perturbations of each reference state. The shifts in free energy obtained by comparison of adjacent simulations and their SIMULATION OF CAVITY FORMATION I 01 I I I I 0-0.06 0.08 0.10 0.12 0.14 FIG.2.-Theoretical free energy of cavity formation (drawn curve) and points derived by perturbation of the Ych = 0.1 nm simulation illustrating inadequacy of extension from simulation to lower values of the cavity radius.accumulated values are given in table 2. We see that connection between Yth = 0.238 and 0.299 involves a rather large shift that may require an intermediate simu- lation. From the present analysis we cannot construct the AH curves with acceptable accuracy. Also the volume allows no definite conclusions; compared to the volume of the cavity itself the indications are that a volume contraction occurs that com- pensates the cavity volume itself. TABLE 2.-ADJUSTMENT OF FREE-ENERGY DIFFERENCES BY PERTURBATIONS OF MD SIMULATIONS connection shifts accumulated between rlh/nm /kJ mol-' shift/kJ rno1-I ~~ ~ ~~ 0.100+0.178 3.5 f0.5 3.5 f0.5 0.178 +0.238 7.3 0.6 10.8 f0.7 0.238 -+0.299 0.299 -f 0.317 11.1 -+0.6 4.5 -+ 0.6 21.8 f0.9 26.3 f1.1 4.2.STRUCTURAL PROPERTIES The distribution of water molecules around a cavity is characterized by the cavity- water pair distribution gcw(r),and by the water-water pair distribution function on the cavity surface. Water orientations are characterized by orientational distribution functions of OH HH and dipole vectors. Fig. 3 gives the cavity-water pair distribution function for all five cavities studied. It is evident that the cavity of 0.178 nm thermal radius has the most pronounced radial shell structure with a high first-neighbour peak and a well developed second-neighbour peak. Also the 0.299 nm cavity has relatively better resolved structure. Considering the orientational distributions for one simulation (rth = 0.299 nm) we J.P. M. POSTMA H. J. C. BERENDSEN AND J. R. HAAK observe a higher occurrence of OH orientation of both 0" (radially outward) and 120° while the dipole is mostly oriented around 70-80" (fig. 4). Very roughly a combin- ation of configurations as given in fig. 5 would account qualitatively for the orien- t ational distribution curves. The pair distribution curve of molecules in the first shell can best be expressed in terms of the distribution of the angle 8 between two molecules of any pair in the shell.4 Fig. 6 shows the distribution of cos 8 in the first solvation shell (rth = 0.299 nm) n fi R I x Y .-.I > 3 c5 rlnm FIG.3.-Cavity-oxygen radial distribution curves for five values of the thermal radius of the cavity (A) 0.100 (B) 0.178 (C) 0.238 (D) 0.299 (E) 0.317 nm.Curves are obtained by spline smoothing with deviation from the data not exceeding 0.05. No smoothing was applied for distances below the first maximum where the shell has been divided into two subshells. The subshell directly on the surface of the cavity shows a pronounced first-neighbour peak and a resolved second- neighbour peak. The first peak corresponds to hydrogen-bonded distances of neighbours. Solvation shell structures are different for different cavity sizes; a more complete analysis of these structures is in progress. 4.3. DYNAMIC PROPERTIES As far as dynamic properties of cavity hydration are concerned we limit our analysis here to the residence-time distribution of molecules in the solvation shells and defer the analysis of reorientational motions to a later publication.The residence- time distribution is a measure of the stability of a solvation shell; it is determined both by the size of the shell and the diffusion constant essentially in radial direction. The residence-time distribution was determined by monitoring what fraction of the molecules present in a shell at t = 0 were present at time t. This was averaged over initial times and corrected for the fraction of molecules that are expected to occur in the shell on the average. This distribution is equivalent to the correlation function of a variable that is equal to 1 when a molecule is in the shell and zero when it is not in the SIMULATION OF CAVITY FORMATION 0.8 0.6 0.4 0.2 c -1 -0.5 0 0.5 1 P!' ! ' ! ' I 114 180150 120 90 60 30 0 inner scale cos 0 outer scale O/" FIG.4.-Probability density of orientation of OH dipole and HH direction with respect to radius from cavity centre to oxygen expressed as distribution over cos 8.Data apply to molecules in first shell (Y < 0.475 nm) from simulation with T(h = 0.299 nm. shell. The resulting curves for the first second and third shell of the 0.299 nm cavity are given in fig. 7. After an initial fast partial decay the distribution decays very slowly reaching half-lives of several picoseconds. The interest in these curves lies in the prediction of an effective diffusion constant FIG.5.-TWO possible orientations of water consistent with the orientational distributions of fig.4. J. P. M. POSTMA H. J. C. BERENDSEN AND J. R. HAAK for the molecules in a shell. The exact calculation of the residence time distribution in a spherical shell with a non-uniform distribution is prohibitive but a good approxi- mation is possible. For this we consider a planar slab with thickness a uniformly filled with particles in an infinite medium. The particles diffuse with a diffusion V!'! I ! I 1 1 I4 180150 120 90 60 30 0 inner scale cos 8 outer scale S/" FIG.6.-Pair correlation of water molecules on cavity surface for rth = 0.299 nm. G(cos 0)dcos 8 gives the fraction of the total number of neighbours in the indicated shell that are in an angular range (cos 8 cos 8 + dcos S) given the presence of one molecule at 8 = 0".(-) rcav-ox < 0.4 nm; (--) 0.4 < rlnm < 0.475.constant D through the slab and the medium. The fraction P(t) of the particles originally present in the slab that are still in the slab at time t is given by P(t) = 2adnDt fixExf exp (w2). In the case that one of the walls of the slab is reflecting the lower bound of the second integral should be -a. Solving the integrals we obtain P(z) = 1 -d(z/n)[I -J(z-+)] where J(x) is the normalised integral of the error function also indicated as which is available in tabulated form.23 The variable z is given by z = Dt/a2 for a reflecting slab and z = 4Dtla2 for an open slab. SIMULATION OF CAVITY FORMATION 0 1 2 3 4 5 t ime/ps FIG.7.-Residence-time distributions in (a)first (r< 0.475 nm) (b)second (0.475<r/nm<0.76) and (c)third (0.76<r/nm<0.9) shell of cavity with thermal radius of 0.299 nm.Drawn curves give the fractionP(t) of molecules present in the shell at t = 0 that are also present at time t. P(t)is corrected for the probability of findinga molecule in the shell at t =co due to the finite size of the system. Broken curves are theoretical residence time distributions based on a diffusion model. Applying eqn (23a) and (23c) for the first solvation shell and eqn (23b) and (23d) for the second and third solvation shells we can determine the diffusion constant. The thickness a was chosen as 0.175 0.285 and 0.140 nm for the three shells.Using values for D of 8 x 6 x lov5and 5 x cm2s-' for the first second and third shell respectively we obtain the broken curves given in fig. 7. The measured curve agrees very well with the theoretical one. The values of the diffusion constants in the three shells do not differ significantly from the bulk value of (6.8 0.7) x cm2 s-' of the simulation. The conclusion must be that at least the translational dynamics are not slower in the solvation shell than in the bulk liquid. 5. DISCUSSION We limit our discussion to the comparison with scaled-particle theory; a discus-sion of structural and dynamic properties requires a more complete analysis of the data. We have shown that it is quite possible to obtain values for the free energy of cavity formation by scaling a cavity up from zero to its size through a limited number of discrete simulations.It is much more difficult to obtain enthalpy and entropy values with comparable accuracy. It is of special interest to compare our " experimental " AGc values with those predicted from scaled-particle theory. There are two difficulties (a) what is the hard-sphere radius of a cavity with given thermal radius and (b)what is the hard- 67 J. P. M. POSTMA H. J. C. BERENDSEN AND J. R. HAAK sphere radius a of the water molecule? To the first question there is a unique answer the theoretical exact relation for small hard cavities in a hard-sphere solvent [eqn (3)] is equivalent to the exact relation for soft cavities [eqn (lo)] when the hard-sphere radius rhs is related to the thermal radius as rhs = 1.0704 rth.(24) The numerical factor is equal to [l-(3/4)]1’3. This relation we shall apply throughout. The hard-sphere radius a of a water molecule cannot be determined independently. When a = 0.2875 nm is taken using the density of our water model (0.983 g ~m-~) SPT gives the values of AG as shown in fig. 1. The agreement with “experiment ” is excellent. This value of a is slightly (but significantly) larger than the values of 0.275 or 0.277 nm used in the literature. The parameter y has a value of 0.4092. For the coefficients K in eqn (4) neglecting the term in r3 we find KO= 6.651 (254 Kl = -110.39 (25b) K2 = 507.82. (254 The excellent agreement supports the applicability of scaled-particle theory to water provided a reliable estimate for the effective hard-sphere radius can be made.The quadratic term in the expansion of AG [eqn (4)] relates to the surface tension y. Thus we can derive the surface tension from our simulations. We find y = K2/(4nNAv)kJ nm-2 (26) or y = 0.067 N m-’ which agrees closely with the experimental value at 305 K of 0.071 N m-l. We gratefully acknowledge the generous support of the Computer Centre of the University of Groningen. This research was supported by the Foundation for Chemical Research (S.O.N.) under the auspices of the Netherlands Organisation for the Advancement of Pure Research (Z.W.O.). See for a review H. Reiss Adv. Chem. Phys. 1965 lX,1 and references quoted therein. ’R. A. Pierotti J. Phys.Chem. 1965 69 281. R. A. Pierotti Chem. Rev. 1976 76 717. F. H. Stillinger J. Solution Chem. 1973 2 141. H. M. Neumann J. Solution Chem. 1977 6 33. H. DeVoe J. Am. Chem. SOC.,1976,98 1724. ’M. Lucas J. Phys. Chem. 1976,80 359. M. Lucas and R. Bury J. Phys. Chem. 1976 80 999. N. Morel-Desrosiers and J-P. Morel Can. J. Chem. 1981 59 1. lo V. G. Dashevsky and G. N. Sarkisov Mol. Phys. 1974,27 1271. A. Geiger A. Rahman and F. H. Stillinger J. Chem. Phys. 1979 70,263. l2 C. Pangali M. Rao and B. J. Berne J. Chem. Phys. 1979 71 2975; 2982. l3 G. Alagona and A. Tani J. Chern. Phys. 1980,72 580. l4 K. Nakanishi S. Okazaki K. Ikari and H. Touhara Chem. Phys. Lett. 1981 84 428. l5 H. J. C. Berendsen J. P. M. Postma W. F. van Gunsteren and J. Hermans in Intermolecular Forces ed.B. Pullman (D. Reidel Dordrecht 1981) pp. 331-342. l6 J. 0. Hirschfelder C. F. Curtiss and R. B. Bird Molecular Theory of Gases and Liquids (J. Wiley New York 1954) p. 157. D. A. McQuarrie Statistical Mechanics (Harper and Row New York 1976). J. P. M. Postma unpublished results. l9 W. L. Jorgensen J. Am. Chem. SOC. 1981 103 335. ’O J. P. Ryckaert G. Cicotti and H. J. C. Berendsen J. Comput. Phys. 1977,23 327. W. F. van Gunsteren and H. J. C. Berendsen Mol. Phys. 1977,34,1311. ”H. J. C. Berendsen J. P. M. Postma A. Di Nola W. F. van Gunsteren and J. R. Haak to be published. 23 M. Abramowitz and I. A. Stegun Handbook of MathematicalFunctions (Dover Publication Inc. New York 1965).