- Email: [email protected]
A molecular dynamics study of the TIP4P model of water
Volume 121. number 3
A MOLECULAR M
8 November 1985
CHEMICAL PHYSICS LETTERS
DYNAMICS
STUDY
OF THE
TIP4P
MODEL
OF WATER
FERR%RIO
and A
TANI
Rtzrened 26 June 1985. m final form 16 Augusl
1985
1_ Introduction In the last decade a large number of pairwise addrtlve models of the mtermolecular potential has been developed for water [I] III an attempt to pro-de a detalled and reahstic descnption of the hquld. Among them, the TIP4P model, a member of a class of transferable intermolecular potentrals (TIP) for orgamc and inorgamc molecules, was proposed by Jorgensen et al [2] The TIP4P potential for water successfully reproduces the energy and density of liquid water at 1 atm pressure and 25’C but overestimates the height of the first peak of the oxygen+xygen radial distributlon function. The latter feature might be due to the neglect of many-body effects m an effective par mteractlon model that gives correct values for a number of thermodynamic properties. Its dynamical propertles, however, and its behaviour as a function of temperature, are still to be explored_ In view of the mcreasmg use of this potential model m nmulations of aqueous solutions [3,4], there IS a need for a more accurate and extensive investigatron of Its properties_ With this m mind, we are carrymg out a senes of molecular dynamics (MD) simulations, some of them m the Isothermal-isobaric ensemble Given the long 182
range nature of the water-water potential and our Interest m dielectric and orientatronal propertres, a farrly large system (343 molecules) has been adopted This system size has been combined with Ewald summatlon for the electrostatic interactrons and the maxrmum cutoff radius ava.dable_ The latter factor has been shown [5] IO affect crucially the orientatronal distribution of water dipoles In this Letter we report some results from a series of runs at three different temperatures Simulation details are given in sectron 2 and results are discussed in section 3
2 Computational
details
Molecular dynamics calculations have been carned out on a 343 TIP4P-water molecule system m a cubic box with penodic boundary conditions Oxygenoxygen 12-6 mteractlons have been evaluated with a spherical cutoff radius R, equal to half the MD box length L (R, = 10.93 A at T= 25OC) with usual mnunum image cntena Long-range corrections have been considered for both the energy and the virlal sum by replacing the radral drstnbutron functions 0 009-2614/85/S (North-Holland
03.30 0 Elsevier Physics Publishmg
Science Dtision)
Publishers
B.V
Volume 121, number 3
CHEMICAL PHYSICS LETTERS
coefficients that are large at all temperatures The drelectric constant was estrmeted from the average value of the IGrkwood gK factor, according to the relation
(g(R)) with 1 beyond the cutoff distance R, The electrostatrc long-ranged forces and energies have been calculated by the Ewald method [6] _ The drmenslordess convergence parameter (Yhas been optimized with respect to the real space and the reclprocal space cutoff distances We used LY= 5 8 and a real space cutoff equal to half the MD box length, i e the same as for the 12-6 potential. The reciprocal space cutoff Kc was fared by the relatron K& = 27 U2 (Kc = 0.51 A-l at T= 25°C). The equations of motion have been integrated by the generahzed method of constraints [7] _ The trmesteps m the three reported runs were At = 1 5) 2 0, 2.5 fs for the temperatures 75, 25, -30°C, respectively. The combined use of Ewald sums and the method of constramts gave very good energy conservation with essentially no dnft and fluctuations smaller than IO-3%. Equihbration runs of about 2000 timesteps were partly performed at constant temperature with the extended method of NoG [8,9] momtoring the stabllizatlon of the structure and of some thermodynamic properties Reported results were calculated from configuratrons out of a standard rmcrocanomcal traLengths of these equiliJectory, after equtibratlon brium runs are reported in table 1 together with ngnifIcant thermodynamic properties calculated as run averages Experimental values correspond to temperatures gven in parentheses in table 1 and denstres were chosen equal to the experimental values, except for the 2S°C run where Jorgensen’s onginal result was preferred for the sake of corppanson All results are m good agreement Hrlth expenment apart from diffusion Table 1 Thermodynarmc
steps (p/cm3)
3 Results and drscusion The three atomic radial dlstrrbutlon functions, gxy(R), are collected in fig 1 For comparison the Monte Carlo results at T = 25°C of Jorgensen et al [ 151 are shown in the same figure. The agreement between the present MD, at 27” C, and the Monte Carlo curve 1s quantitative, despite differences in system size (343 versus 216 molecules) and m the treatment of long-range electrostatic interactions (Ewald sum versus spherical cutoff). This confirms the well known insensitivity of these functions to these factors [ 5]
300 0 37 9.85 24 6 45
9 92 17.6 a) 2 4 C)
2100 0.978 351 103 9 20
_
(0.997) (298) 0 01
72 2-l
78”) 25 d)
-
8100
10 74 (10 4) OS
A
(28) SO [ll]
Expt
2500
-E (kcal/mol) Cy(cal/mol°C) 105D (cm2/s)
b, Ref
Run 3
Run 2
c) Ref
0.983
0 985
0 985
(PS)
Eupt.
Expt
242 0 06
[lo]
,
where y = 4ap2,/9X-, T, which IS gwen m ref [ 141 for the Ewald sum case The reported specific heats at constant volume are augmented by the estimated vIbrational contributions The values of both dielectric constant and specific heat for run 1 and run 3, still in progress, are ,gven in parentheses, the length of these runs bemg too short to give reliable values for these properties related to the variance rather than to the mean of the relevant distribution. Most calculations were performed on a Could-Se1 32187, part on an IEIM 3081 and part on a CRAY XMP. Relative computer speeds were roughly 4 for IBM FORTRAN VS and 40 for CRAY FORTRAN with respect to the Gould
Run 1
T(K) P &bar)
a) Ref
E - 1 = SygK(Ewald)
properties
number of densty
8 November 1985
(243) 0 001 10.76 17 s=) 02a) - 104 a) [12]
‘1 Ref. [
(10 0) 7.6 (44) 1.1
0 978 (348) 0.001 9 21 16.7 b) 6 0 c) 62b) I -05 d)
131 183
Volume 121. number 3
CHEMICAL PHYSICS LETTERS
ITT-m
1 5.. 1 -.
-5..
2
4
6 rlH
0
Fg_ 1. Radial dlstnbutlon functions for atormc pans Open
circles ref [15],T=298K_Dottedllne Solidhncs
runl.T=242K
run2,T=300K,andruo3,P=351K.
Temperature changes produce mamly quantltatrve effects m the range of distances relevant to nearestnelghbour molecules For example, the trend toward more precise tetrahedral arrangement when temperature decreases can be monitored through the values of coordmation numbers for the first peaks ofgoo(R) (5367 at 351 K), 4515 at 300K and 4-121 at 242K) and ofgoH(R) (0.898,0.928 and 0.970 respectively). The supercooled hqurd, however, shows evidence 184
8 November 1985
of enhanced order also in those peaks that relate to mteractions of second and higher order nearest neighbours This evidence is both quantitative, the larger amphtudes of the osctiatlon ofg&R), and qu&taWe, m particular the oscillations in the region 4.5 to 7 0 A ofgHH(R) and goH(R) hi the latter function the oscillations have opposte phase at 242 and 35 1 K The postion of the maxuna ofgo, ISmdependent of temperature, in contrast to other models of water such as ST2 [ 161, where the second peak slufts to larger R as temperature increases, and MCY [ 171 that shows an opposite trend It should be noted that the calculations in ref [ 171 were performed at a constant denslty of 18.0 cm3/mole so that the Increase of temperature leads to an increase of pressure, with a concurrent shift of the second peak to shorter drstances In fact a remarkable mward shaft of this peak can be observed companng Owicki and Scheraga’s [ 181 results (maximum at 4-75 A) at the constant pressure of 1 atm W&I those [17] obtained at the expzrunental density and hence higher pressure (mazemumat43A) The dlstributlon of par interactron energies (fig. 2), has been calculated for molecular pars whose oxygen-xygen separation does not exceed 3 5 A, m order to build a data base for a detaled analysis of static and dynamic propertres of hydrogen-bonded molecules [19] Th.~schoice of 3 5 .& as distance threshold, although arbitrary, has allowed us to extract from the total distribution all the significant contnbutions, both attractive and repulsive- There is quanhtatlve agreement between Jorgensen’s [15] and our results up to -3 0 and above +3 0 kcal/mole. The overall shape of the distributions suggests the existence of a number of discrete populations, presumably related to hydrogen-bonded and non-bonded pairs, whose relative weight depends on temperature. In this respect, these results differ from those obtained for MCY water at T = 2S°C and P= 1 atm [ 181, which show a continuous transitlon from attractive to repulsive mteractlons. Comparison with other models 1s obscured by the presence of the overwhelming central peak. as this dxtrlbutlon has been calculated for all pairs, irrespective of their separation_ The dependence on tevperature, however, is quahtatlvely the same as that reported for BNS and ST2 [ 161 water. In particular, here too we fmd a nearly mvanant value of interaction energy located at -3.9
Volume
121. number 3
CHEMICAL NEAREST .4
NEIGHBOUR
PHYSICS LETTERS
ENERGY
DISTRl9LlTION
1
t
-45 FG. 2 Pair energy dlsmbution ior Open arcles ref [ 151, T = 298 K
-25
-0.5
15
kcal Imol
newest ncighbour molecules. (1) Run 1, T = 242 K, (2) run 2.T=300K;(3)run3.T=351K
4 0 for ST2 and -3.7 k&/mole for BNS Thrs behaviour of the pau energy drstnbution as a functron of temperature has been mterpretated as an mdrcatron of a “hydrogen-bond rupture mechamsm with a dynamic equrlrbnum between bonded and non-bonded pars” [ 161, i e paxs with interaction energy below and above -3 9 kcal/ mole, in our case. The same mechanism is considered to account for the rsosbestrc point experimentally observed for OH and OD mtramolecular stretchmg bands [20] wrth Raman and infrared spectra at various temperatures_ The reorientational motion of water m the liquid phase can be studied through a family of autocorrelatron functrons (ACF) of the mdrvrdual dipole drrections, versus
Iegendre polynonnal of order k and vector_ We have computed the first fig. 3 shows only the results fork = 1, respect to a coordinate frame fved on the molecule Values of the relaxation time T1, as derived from the slope of the curves in fig. 3, are reported rn table 1.
where Pk IS the H is the dipole four C,(t), but calculated with
8 November 1985
Quahtatrvely, these functions exhrbit the same behavrour observed for BNS or MCY water, i.e. a rebound at short times followed by a diffusive exponential decay. Thus is consrdered to indicate the existence of cages of four tetrahcdrally coordinated molecules that host the tagged water molecule_ The latter 0 PZ_CIUXCZPR?_<~IS~ ‘L-r-,Cr
?
0
1
2
Fig_ 3 Logarithm plot of the dipole nutomrrelatmn (l)Runl,T=242K;(2)run2.T=300K;(;)~n3,T= 351 K
3
upr
function
185
Volume 121, number 3
CHEMICAL PHYSICS LETTERS
undergoes a hindered oscrllation or hbration that marufests Itself in the shortest time part-of C,(r), whtle the slower breakup of the cage 1s associated with the exponentral decay at larger times. Tlus mechanism seems to be conststent w-nh the observed temperature dependence of Cl(f).
Admowledgement We wish to thank W L Jorgensen for providing us with delalled resultsofhis calculation and the Istituto di Chimica Quantistica, CNR, Plsa, for a generous allocatron of computer time on its Gould 32/87
References [l]
J R. Reimers. R-0. Watts and hi L Klcm. Chem
Phys
64 (1982) 9.5, and references therein. [Z] W L Jorgensen. 1 Chandrasekhar, J D. Madura, R W. lmpey and M L Klem, J. Chcm Phys 79 (1983) 926; W L Jogcnscn, J Am Chcm Sot 103 (1981) 335 [3] W-L. Jorgensen and CJ. Swenson, J Am Chem. Sot. 107 (1985) 1489. and references theran [4] G Alagona, C Ghlo and P Kollman. J Am Chem Sot. (1985). to be pubhshed
186
8 November 1985
[S] T A Andrea, W C Swope and H C. Andcrsen. J Chem Phys 79 (1983) 4576 (61 E R Cowley. G Iacucci, M L Klein and I R.McDonaJd. Phys Rev. B14 (1976) 1758 [7] G Clccotti, JP Ryckaert and M. rerrano, Mol Phys 47 (1982) 1253. [8] S NOSE, J. Cbem. Phys_ 81 (1984) 511. [9] M Ferrano and JP_ Ryckaert. Mol. Phys 54 (198.5) 587 [lo] C A AngeJl. ~1: Water. a cornprehensIve treatise. Vol 7, ed. F Franks (Plenum Press, New York, 1982) ch. 1 [ 111 J B. Hasted. m- Water. a comprehensive treatise, Vol 7, ed F Franks (Plenum Press, New York. 1982) ch 7. [12] K. Krynicki, CD Green and D W. Sawyer, Faraday Drscussions Chem. Sot. 66 (1978) 199 [ 131 H.G Hertz and M.D. Zeidler, m The hydrogen bond, Vol 3. eds. P. Schuster, G Zundcl and C Sandorfy (North-Holland, Amsterdam, 1976) ch 21. [14] G Stell, G_N Patey and J S Hoyc. Advan Chem. Phys 48 (1981) 183_ [ 151 W L Jorgensen. pnvate commumcatron [16] F H. Stinger and A Rahmar1.J Chem Phys 61 (1974) 4973 [17] R W Impey. M L Klnn and I R McDonald, J Chem Phys 74 (1981) 647 [I81 J.C. Owicki and H-A. Scheraga. J- Am. Chem Sot 99 (1977) 7403. [ 191 M rerrano and A Tam, work m progress [20] GE Walrafen. m: Water, a comprchcnsive treatise, Vol 1. ed P Frarks (Plenum Press,New York, 1972) ch 5.
A MOLECULAR M
8 November 1985
CHEMICAL PHYSICS LETTERS
DYNAMICS
STUDY
OF THE
TIP4P
MODEL
OF WATER
FERR%RIO
and A
TANI
Rtzrened 26 June 1985. m final form 16 Augusl
1985
1_ Introduction In the last decade a large number of pairwise addrtlve models of the mtermolecular potential has been developed for water [I] III an attempt to pro-de a detalled and reahstic descnption of the hquld. Among them, the TIP4P model, a member of a class of transferable intermolecular potentrals (TIP) for orgamc and inorgamc molecules, was proposed by Jorgensen et al [2] The TIP4P potential for water successfully reproduces the energy and density of liquid water at 1 atm pressure and 25’C but overestimates the height of the first peak of the oxygen+xygen radial distributlon function. The latter feature might be due to the neglect of many-body effects m an effective par mteractlon model that gives correct values for a number of thermodynamic properties. Its dynamical propertles, however, and its behaviour as a function of temperature, are still to be explored_ In view of the mcreasmg use of this potential model m nmulations of aqueous solutions [3,4], there IS a need for a more accurate and extensive investigatron of Its properties_ With this m mind, we are carrymg out a senes of molecular dynamics (MD) simulations, some of them m the Isothermal-isobaric ensemble Given the long 182
range nature of the water-water potential and our Interest m dielectric and orientatronal propertres, a farrly large system (343 molecules) has been adopted This system size has been combined with Ewald summatlon for the electrostatic interactrons and the maxrmum cutoff radius ava.dable_ The latter factor has been shown [5] IO affect crucially the orientatronal distribution of water dipoles In this Letter we report some results from a series of runs at three different temperatures Simulation details are given in sectron 2 and results are discussed in section 3
2 Computational
details
Molecular dynamics calculations have been carned out on a 343 TIP4P-water molecule system m a cubic box with penodic boundary conditions Oxygenoxygen 12-6 mteractlons have been evaluated with a spherical cutoff radius R, equal to half the MD box length L (R, = 10.93 A at T= 25OC) with usual mnunum image cntena Long-range corrections have been considered for both the energy and the virlal sum by replacing the radral drstnbutron functions 0 009-2614/85/S (North-Holland
03.30 0 Elsevier Physics Publishmg
Science Dtision)
Publishers
B.V
Volume 121, number 3
CHEMICAL PHYSICS LETTERS
coefficients that are large at all temperatures The drelectric constant was estrmeted from the average value of the IGrkwood gK factor, according to the relation
(g(R)) with 1 beyond the cutoff distance R, The electrostatrc long-ranged forces and energies have been calculated by the Ewald method [6] _ The drmenslordess convergence parameter (Yhas been optimized with respect to the real space and the reclprocal space cutoff distances We used LY= 5 8 and a real space cutoff equal to half the MD box length, i e the same as for the 12-6 potential. The reciprocal space cutoff Kc was fared by the relatron K& = 27 U2 (Kc = 0.51 A-l at T= 25°C). The equations of motion have been integrated by the generahzed method of constraints [7] _ The trmesteps m the three reported runs were At = 1 5) 2 0, 2.5 fs for the temperatures 75, 25, -30°C, respectively. The combined use of Ewald sums and the method of constramts gave very good energy conservation with essentially no dnft and fluctuations smaller than IO-3%. Equihbration runs of about 2000 timesteps were partly performed at constant temperature with the extended method of NoG [8,9] momtoring the stabllizatlon of the structure and of some thermodynamic properties Reported results were calculated from configuratrons out of a standard rmcrocanomcal traLengths of these equiliJectory, after equtibratlon brium runs are reported in table 1 together with ngnifIcant thermodynamic properties calculated as run averages Experimental values correspond to temperatures gven in parentheses in table 1 and denstres were chosen equal to the experimental values, except for the 2S°C run where Jorgensen’s onginal result was preferred for the sake of corppanson All results are m good agreement Hrlth expenment apart from diffusion Table 1 Thermodynarmc
steps (p/cm3)
3 Results and drscusion The three atomic radial dlstrrbutlon functions, gxy(R), are collected in fig 1 For comparison the Monte Carlo results at T = 25°C of Jorgensen et al [ 151 are shown in the same figure. The agreement between the present MD, at 27” C, and the Monte Carlo curve 1s quantitative, despite differences in system size (343 versus 216 molecules) and m the treatment of long-range electrostatic interactions (Ewald sum versus spherical cutoff). This confirms the well known insensitivity of these functions to these factors [ 5]
300 0 37 9.85 24 6 45
9 92 17.6 a) 2 4 C)
2100 0.978 351 103 9 20
_
(0.997) (298) 0 01
72 2-l
78”) 25 d)
-
8100
10 74 (10 4) OS
A
(28) SO [ll]
Expt
2500
-E (kcal/mol) Cy(cal/mol°C) 105D (cm2/s)
b, Ref
Run 3
Run 2
c) Ref
0.983
0 985
0 985
(PS)
Eupt.
Expt
242 0 06
[lo]
,
where y = 4ap2,/9X-, T, which IS gwen m ref [ 141 for the Ewald sum case The reported specific heats at constant volume are augmented by the estimated vIbrational contributions The values of both dielectric constant and specific heat for run 1 and run 3, still in progress, are ,gven in parentheses, the length of these runs bemg too short to give reliable values for these properties related to the variance rather than to the mean of the relevant distribution. Most calculations were performed on a Could-Se1 32187, part on an IEIM 3081 and part on a CRAY XMP. Relative computer speeds were roughly 4 for IBM FORTRAN VS and 40 for CRAY FORTRAN with respect to the Gould
Run 1
T(K) P &bar)
a) Ref
E - 1 = SygK(Ewald)
properties
number of densty
8 November 1985
(243) 0 001 10.76 17 s=) 02a) - 104 a) [12]
‘1 Ref. [
(10 0) 7.6 (44) 1.1
0 978 (348) 0.001 9 21 16.7 b) 6 0 c) 62b) I -05 d)
131 183
Volume 121. number 3
CHEMICAL PHYSICS LETTERS
ITT-m
1 5.. 1 -.
-5..
2
4
6 rlH
0
Fg_ 1. Radial dlstnbutlon functions for atormc pans Open
circles ref [15],T=298K_Dottedllne Solidhncs
runl.T=242K
run2,T=300K,andruo3,P=351K.
Temperature changes produce mamly quantltatrve effects m the range of distances relevant to nearestnelghbour molecules For example, the trend toward more precise tetrahedral arrangement when temperature decreases can be monitored through the values of coordmation numbers for the first peaks ofgoo(R) (5367 at 351 K), 4515 at 300K and 4-121 at 242K) and ofgoH(R) (0.898,0.928 and 0.970 respectively). The supercooled hqurd, however, shows evidence 184
8 November 1985
of enhanced order also in those peaks that relate to mteractions of second and higher order nearest neighbours This evidence is both quantitative, the larger amphtudes of the osctiatlon ofg&R), and qu&taWe, m particular the oscillations in the region 4.5 to 7 0 A ofgHH(R) and goH(R) hi the latter function the oscillations have opposte phase at 242 and 35 1 K The postion of the maxuna ofgo, ISmdependent of temperature, in contrast to other models of water such as ST2 [ 161, where the second peak slufts to larger R as temperature increases, and MCY [ 171 that shows an opposite trend It should be noted that the calculations in ref [ 171 were performed at a constant denslty of 18.0 cm3/mole so that the Increase of temperature leads to an increase of pressure, with a concurrent shift of the second peak to shorter drstances In fact a remarkable mward shaft of this peak can be observed companng Owicki and Scheraga’s [ 181 results (maximum at 4-75 A) at the constant pressure of 1 atm W&I those [17] obtained at the expzrunental density and hence higher pressure (mazemumat43A) The dlstributlon of par interactron energies (fig. 2), has been calculated for molecular pars whose oxygen-xygen separation does not exceed 3 5 A, m order to build a data base for a detaled analysis of static and dynamic propertres of hydrogen-bonded molecules [19] Th.~schoice of 3 5 .& as distance threshold, although arbitrary, has allowed us to extract from the total distribution all the significant contnbutions, both attractive and repulsive- There is quanhtatlve agreement between Jorgensen’s [15] and our results up to -3 0 and above +3 0 kcal/mole. The overall shape of the distributions suggests the existence of a number of discrete populations, presumably related to hydrogen-bonded and non-bonded pairs, whose relative weight depends on temperature. In this respect, these results differ from those obtained for MCY water at T = 2S°C and P= 1 atm [ 181, which show a continuous transitlon from attractive to repulsive mteractlons. Comparison with other models 1s obscured by the presence of the overwhelming central peak. as this dxtrlbutlon has been calculated for all pairs, irrespective of their separation_ The dependence on tevperature, however, is quahtatlvely the same as that reported for BNS and ST2 [ 161 water. In particular, here too we fmd a nearly mvanant value of interaction energy located at -3.9
Volume
121. number 3
CHEMICAL NEAREST .4
NEIGHBOUR
PHYSICS LETTERS
ENERGY
DISTRl9LlTION
1
t
-45 FG. 2 Pair energy dlsmbution ior Open arcles ref [ 151, T = 298 K
-25
-0.5
15
kcal Imol
newest ncighbour molecules. (1) Run 1, T = 242 K, (2) run 2.T=300K;(3)run3.T=351K
4 0 for ST2 and -3.7 k&/mole for BNS Thrs behaviour of the pau energy drstnbution as a functron of temperature has been mterpretated as an mdrcatron of a “hydrogen-bond rupture mechamsm with a dynamic equrlrbnum between bonded and non-bonded pars” [ 161, i e paxs with interaction energy below and above -3 9 kcal/ mole, in our case. The same mechanism is considered to account for the rsosbestrc point experimentally observed for OH and OD mtramolecular stretchmg bands [20] wrth Raman and infrared spectra at various temperatures_ The reorientational motion of water m the liquid phase can be studied through a family of autocorrelatron functrons (ACF) of the mdrvrdual dipole drrections, versus
Iegendre polynonnal of order k and vector_ We have computed the first fig. 3 shows only the results fork = 1, respect to a coordinate frame fved on the molecule Values of the relaxation time T1, as derived from the slope of the curves in fig. 3, are reported rn table 1.
where Pk IS the H is the dipole four C,(t), but calculated with
8 November 1985
Quahtatrvely, these functions exhrbit the same behavrour observed for BNS or MCY water, i.e. a rebound at short times followed by a diffusive exponential decay. Thus is consrdered to indicate the existence of cages of four tetrahcdrally coordinated molecules that host the tagged water molecule_ The latter 0 PZ_CIUXCZPR?_<~IS~ ‘L-r-,Cr
?
0
1
2
Fig_ 3 Logarithm plot of the dipole nutomrrelatmn (l)Runl,T=242K;(2)run2.T=300K;(;)~n3,T= 351 K
3
upr
function
185
Volume 121, number 3
CHEMICAL PHYSICS LETTERS
undergoes a hindered oscrllation or hbration that marufests Itself in the shortest time part-of C,(r), whtle the slower breakup of the cage 1s associated with the exponentral decay at larger times. Tlus mechanism seems to be conststent w-nh the observed temperature dependence of Cl(f).
Admowledgement We wish to thank W L Jorgensen for providing us with delalled resultsofhis calculation and the Istituto di Chimica Quantistica, CNR, Plsa, for a generous allocatron of computer time on its Gould 32/87
References [l]
J R. Reimers. R-0. Watts and hi L Klcm. Chem
Phys
64 (1982) 9.5, and references therein. [Z] W L Jorgensen. 1 Chandrasekhar, J D. Madura, R W. lmpey and M L Klem, J. Chcm Phys 79 (1983) 926; W L Jogcnscn, J Am Chcm Sot 103 (1981) 335 [3] W-L. Jorgensen and CJ. Swenson, J Am Chem. Sot. 107 (1985) 1489. and references theran [4] G Alagona, C Ghlo and P Kollman. J Am Chem Sot. (1985). to be pubhshed
186
8 November 1985
[S] T A Andrea, W C Swope and H C. Andcrsen. J Chem Phys 79 (1983) 4576 (61 E R Cowley. G Iacucci, M L Klein and I R.McDonaJd. Phys Rev. B14 (1976) 1758 [7] G Clccotti, JP Ryckaert and M. rerrano, Mol Phys 47 (1982) 1253. [8] S NOSE, J. Cbem. Phys_ 81 (1984) 511. [9] M Ferrano and JP_ Ryckaert. Mol. Phys 54 (198.5) 587 [lo] C A AngeJl. ~1: Water. a cornprehensIve treatise. Vol 7, ed. F Franks (Plenum Press, New York, 1982) ch. 1 [ 111 J B. Hasted. m- Water. a comprehensive treatise, Vol 7, ed F Franks (Plenum Press, New York. 1982) ch 7. [12] K. Krynicki, CD Green and D W. Sawyer, Faraday Drscussions Chem. Sot. 66 (1978) 199 [ 131 H.G Hertz and M.D. Zeidler, m The hydrogen bond, Vol 3. eds. P. Schuster, G Zundcl and C Sandorfy (North-Holland, Amsterdam, 1976) ch 21. [14] G Stell, G_N Patey and J S Hoyc. Advan Chem. Phys 48 (1981) 183_ [ 151 W L Jorgensen. pnvate commumcatron [16] F H. Stinger and A Rahmar1.J Chem Phys 61 (1974) 4973 [17] R W Impey. M L Klnn and I R McDonald, J Chem Phys 74 (1981) 647 [I81 J.C. Owicki and H-A. Scheraga. J- Am. Chem Sot 99 (1977) 7403. [ 191 M rerrano and A Tam, work m progress [20] GE Walrafen. m: Water, a comprchcnsive treatise, Vol 1. ed P Frarks (Plenum Press,New York, 1972) ch 5.