- Email: [email protected]
E water at 277 and 300 K
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
3 December 1993
The viscosity of SPC and SPC/E water at 277 and 300 K Paul E. Smith and Wilfred F. van Gunsteren Departmentof PhysicalChemistry SwissFederal Instituteof TechnologyZiirich,ETH-Zentrum, 809.2Zurich, Switzer[and Received 20 September 1993
The shear viscosity of the SPC and SPC/E water models at 277 and 300 K has been determined using equilibrium molecular dynamics. Good agreement with experiment was obtained for the SPC/E model at both temperatures. The viscosity of the SPC model was lower than the experimental value at both temperatures and displayed only a small temperature dependence The effect of using a reaction field treatment for the long-range electrostatic interaction was also investigated and found to be negligible.
1. Introduction Transport coefficients (diffusion constant, bulk and shear viscosity, thermal conductivity) are important properties which characterize the macroscopic properties of molecular systems [ 1,2 1. Any proposed microscopic model for an atomic system should be able to reproduce the corresponding experimentally observed transport coefficients. For instance, in molecular dynamics simulations, one of the criteria for an acceptable water model is that the model correctly reproduces the self-diffusion constant of pure water. Formulae exist for the determination of all the transport properties of a particular model from atomic level simulations [ 11. However, while the selfdiffusion constant is a single-molecule property and is therefore relatively easy to determine, the viscosity and thermal conductivity are properties of the whole system and their determination is correspondingly more difficult. The poor statistics available for these properties means that their determination requires either very long equilibrium simulations, or the use of nonequilibrium molecular dynamics. Experimentally, the shear viscosity of water has been well studied. In contrast, although many models have been proposed for the study of liquid water, which reproduce a number of properties of liquid water, the viscosity of these models has not (to our knowledge) been determined. This is due, in part, to the poor statistics inherent to these methods, and to
the confusion concerning the most efficient method for the determination of shear viscosity from an equilibrium simulation [ 2-5 1. Here we report the determination of the zero-frequency shear viscosity for the SPC [ 61 and SPC/E [ 71 models of water at two different temperatures, 277 and 300 K, using long equilibrium molecular dynamics simulations. The effect of including long-range interactions using the reaction field technique was also investigated.
2. Method A total of eight molecular dynamics simulations were performed. The SPC [6] and SPC/E [ 71 models of water were simulated at 277 and 300 K using a reaction field potential [ 81 with a permittivity of either eRr= 1 (no reaction field) and eRF= cc. All simulations involved a system of 5 12 waters in a cubic periodic box. The systems were initially equilibrated in the NPT ensemble for 50 ps. The average box lengths from these simulations were then used as the box dimensions for a series of NVT simulations. A further 50 ps equilibration was performed in the NVT ensemble. Each system was then simulated for a further 1000 ps. All simulations were performed using the GROMOS program [ 9 1, with a cutoff distance of Rc = 0.9nm. Temperature regulation was achieved with the weak coupling method [ lo] and a coupling constant of 0.1 ps. Bond lengths and angles were constrained using SHARE [ 111 and
0009-261419316 06.00 0 1993 Elsevier Science Publishers B.V. All rights reserved.
315
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
a relative tolerance of 10m5.The time step was 2 fs. The shear viscosity is related to the rate of momentum transport through the system. The shear viscosity can be obtained from simulation by either Green-Kubo- or Einstein-type expressions. The Green-Kubo approach relates the shear viscosity q to the time autocorrelation function of the off-diagonal elements of the pressure tensor [ 1,12 1,
3 December 1993
t
AG(t)=G(t)-G(O)=
j”i;(t’) dt’,
(5)
0
where the time derivative of G(t) is given by .
G(t)=
’
7
=Pa/3(0 , (1)
where Y is the volume of the system, kB is the Boltzmann constant, T is the absolute temperature, t is time, and Pap is the a/3 component of the pressure tensor (cr,+x, y, z; cr#B) given by
PC@(t) = _: x
5 ( i=l
+ mi
“d ‘= f.“, 2k,Tdt
5 icj
Gja(wj/.d~)
3
(2)
>
(
[G(t)-G(0)12) ,
(3)
where G( t ) is given by
ijglria(t)Pi,dt)
and ri( t) and pi( t) are the position and momentum of the centre of mass of molecule i at time t, respectively. The angular brackets in eqs. ( 1) and (3) denote an average over time origins. Unfortunately, G(t) is not continuous under periodic boundary conditions. However, one may write [ 161 316
and& is the total force on atom i. Then one has an Einstein-type relationship which is independent under periodic boundary conditions. This is the approach we have adopted here. During the simulations, each element of the pressure tensor was saved at every step for analysis. We have used these elements of the pressure tensor to calculate the viscosity using eqs. (3), (5) and (6).
Pia(f>Pi#(t>
pr( t) is the centre of mass momentum of molecule i (mass mi) at time t, rij= ri - rj is the intermolecular distance vector, and& is the force on molecule i by molecule j, and the summation is over all N molecules. However, the above correlation function often contains significant long time tails [ 13 1, which have to be integrated correctly to obtain a reliable value of ?/. To avoid the problem of long time tails many transport properties are determined using the equivalent Einstein expression. Helfand derived the following, generally accepted [ 35 1, Einstein relationship for the shear viscosity [ 14,15 ] :
G(t)=
(6)
3. Results The time development of AG(t)’ from one of the simulations is shown in fig. 1. The other simulations produced data of a similar quality. After an initial nonlinear regime of about 5 ps, the curves become linear in time. The shear viscosity was calculated 2000 ,
I
I
J
1500 ?-.
a
“E
s
P 1000 a n “E
s? ”
500
0
0
2
6
4 Time
8
10
(ps)
Fig. 1. The mean-square displacement of G(t) as a function of time obtained by integrating the off-diagonal elements of the pressure tensor. The curves are for the SPC/E model simulated at 300 K without reaction field. The thin lines correspond to the three independent off-diagonal elements of the pressure tensor, and the thick line is the resulting average.
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
from the slope of these curves obtained between 5 10 ps. Errors were determined from the variation between the results for the different off-diagonal components. The final values are presented in table 1. The results for the SPC models are not so good. The viscosity is too low at both temperatures and there is too small a variation of viscosity with temperature. The SPC/E model, however, performs very well with excellent agreement at 300 K and acceptable agreement at 277 K, although the viscosity is slightly too low. The above results reflect nicely the differences between the SPC and SPC/E water models. The SPC/ E model has slightly larger partial atomic charges than the SPC model (qo= -0.8476 versus -0.8200 e). Hence, the SPC/E model has a higher density, a lower diffusion constant and a larger viscosity than the SPC model. For all three properties the SPC/E model gives better agreement with experiment. It is interesting to note that using a reaction field treatment of the long-range interactions did not appear to have any significant effect on the viscosity. However, the reaction field did affect the self-diffusion constant of water as observed in the simulations (see table 1). It is known that there is an in-
and
timate relationship [ 19,201 between viscosity and self-diffusion (Da 1 /q ) , and therefore one might expect a significant variation in viscosity with eRF.Table 1 shows the product Dr7u/k,J, where d is a molecular diameter. This product should be constant if Stokes law is obeyed [ 191. Experimentally, it is constant over the temperature range covered here. However, the simulation values vary by as much as 30% between different models and reaction fields, although this is probably due to the sizeable errors estimated for the calculations. In fig. 2 we display the pressure autocorrelation functions calculated for the SPC/E model at 277 and 300 K without a reaction field treatment. Both the diagonal and off-diagonal components of the pressure tensor are shown in fig. 2a. The diagonal elements display a large initial loss in correlation, followed by a very slow decay over many picoseconds. The off-diagonal elements display an even faster initial decay and a smaller long-time correlation. Both types of correlation function contain some periodicity; a periodicity common to both functions of about 0.10 ps, and a smaller period for the off-diagonal elements of about 0.05 ps. In fig. 2b the difference in the correlation functions with tempera-
Table 1 Molecular dynamics results
T(K)
fIlF
Box length L
Density
Diffusion m*/s)
SPC 277 300
1 co 1 W
SPC/E 277
1 W
300
exp. 0) 277 300
1 co
Viscosity
Pressure
(cp) *)
(W/mol/nm3)
Dqu b)'kBT
P
2.49015 2.50756 2.49975 2.52286
0.991 0.971 0.980 0.953
2.80 3.59 3.94 5.28
0.79( 16) 0.80(4) 0.54(9) 0.58(9)
-3.27 0.85 1.73 0.11
0.223 0.292 0.199 0.289
2.47209 2.49107 2.48133 2.50257
1.013 0.990 1.002 0.976
1.62 1.75 2.68 3.18
1.39(6) 1.40(21) 0.82(9) 0.91(7)
-1.65 -0.28 -0.28 -2.26
0.226 0.247 0.204 0.271
1.000 0.996
1.1 *I 2.2 *)
1.57 0.85
0.174 0.174
a) Standard deviations are given in parentheses. b, The molecular diameter u was obtained from 2 x(IY/~)~= L 3/5 12. Cl Ref. 1171. *‘Ref. 1181.
317
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
3 December 1993
found to be negligible. The results suggest that the SPC model should not be used for low-temperature studies, or for simulations where hydrodynamical effects are critical.
2000
1500
1000
Acknowledgement
500
0
-500 t 0.0
0.2
0.4
0.6
0.6
1.0
1000
We thank Ilario G. Tironi for valuable discussions and gratefully acknowledge the generous financial support of Zeneca Ltd which made this study possible.
600
600
References 400
[ 1] R. Zwanzig, Ann. Rev. Phys. Chem. 16 (1965) 67. [ 21 D.A. McQuarrie, Statistical mechanics (Harper and Row,
200
0
.200' 0.0
a 0.2
s 0.4
' 0.6
s 0.6
' 1.0
Time(ps)
Fig. 2. Pressure time autocorrelation functions for the SPC/E model simulated without reaction field. In (a) the upper curves correspond to the diagonal elements (xx, yy and zz), and the lower curves correspond to the off-diagonal elements (xy, xz and yz) at 300 K. In (b) the thick lines correspond to 300 K and the thin lines to 277 K. The top two curves correspond to the xx component of the pressure tensor and the bottom two curves to the xy component.
ture is displayed. Differences in the correlation functions appear at long times and are greater for the diagonal elements. Integration of the off-diagonal elements to 10 ps gave viscosities which were essentially within the error of those calculated using the Einstein-type relationship. Integration of the diagonal elements (to give the bulk viscosity) was not possible.
4. Conclusions The zero-frequency shear viscosity of SPC and SPC/E water has been determined using long equilibrium molecular dynamics simulations. The SPC/ E model reproduces the experimental viscosity at both 277 and 300 IS. The SPC model does not. The effect of long-range interaction on the viscosity was 318
New York, 1976). [3]A.A.ChialvoandP.G.Debenedetti,Phys.Rev.A43 (1991) 4289. [4] A.A. Chialvo, P.T. Cummings and D.J. Evans, Phys. Rev. E 47 (1993) 1702. [5]M.P.AllenandA.J. Masters,Mol.Phys. 79 (1993) 435. [6] H.J.C. Beret&en, J.P.M. Postma, W.F.vanGunsterenand J. Hermans, in: Intermolecular forces. Interaction models for water in relation to protein hydration, ed. B. Pullman (Reidel, Dordrecht, 1981) pp. 331-342. [7] H.J.C. Berendsen, J.R. GrigeraandT.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [8] J.A. BarkerandR.0. Watts, Mol. Phys. 26 (1973) 789. [9] W.F. van Gunsteren and H.J.C. Berendsen, Groningen molecular simulation (GRO-MOS) Library Manual, Biomos, Groningen (1987). [ lo] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola and J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [ 111 J.-P. Ryckaert, G. Ciccotti and H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. 1121 F.C. Andrews, J. Chem. Phys. 47 (1967) 3161. [ 131 B.J. Alder, in: Molecular dynamics simulation of statistical mechanical systems. Molecular dynamics simulations, eds. G. Ciccotti and W.G. Hoover (North-Holland Amsterdam, 1986) pp. 66-80. [14]E.Helfand,Phys.Rev. 119 (1960) 1. [ 151 D. Gass, J. Chem. Phys. 51 (1969) 4560. [ 16 ] B.J. Alder, D.M. Gass and T.E. Wainwright, J. Chem. Phys. 53 (1970) 3813. [ 171 R.C. Weast, Handbook of chemistry and physics (CRC Press, Boca Raton, 1975). [ 18 ] K.T. Gillen, DC. Douglas and M.J.R. Ho&, J. Chem. Phys. 57 (1972) 5117. [ 191 R. Zwanzig, J. Chem. Phys. 79 (1983) 4507. [20] R. Zwanzig and AK. Harrison, J. Chem. Phys. 83 (1985) 5861.
CHEMICAL PHYSICS LETTERS
3 December 1993
The viscosity of SPC and SPC/E water at 277 and 300 K Paul E. Smith and Wilfred F. van Gunsteren Departmentof PhysicalChemistry SwissFederal Instituteof TechnologyZiirich,ETH-Zentrum, 809.2Zurich, Switzer[and Received 20 September 1993
The shear viscosity of the SPC and SPC/E water models at 277 and 300 K has been determined using equilibrium molecular dynamics. Good agreement with experiment was obtained for the SPC/E model at both temperatures. The viscosity of the SPC model was lower than the experimental value at both temperatures and displayed only a small temperature dependence The effect of using a reaction field treatment for the long-range electrostatic interaction was also investigated and found to be negligible.
1. Introduction Transport coefficients (diffusion constant, bulk and shear viscosity, thermal conductivity) are important properties which characterize the macroscopic properties of molecular systems [ 1,2 1. Any proposed microscopic model for an atomic system should be able to reproduce the corresponding experimentally observed transport coefficients. For instance, in molecular dynamics simulations, one of the criteria for an acceptable water model is that the model correctly reproduces the self-diffusion constant of pure water. Formulae exist for the determination of all the transport properties of a particular model from atomic level simulations [ 11. However, while the selfdiffusion constant is a single-molecule property and is therefore relatively easy to determine, the viscosity and thermal conductivity are properties of the whole system and their determination is correspondingly more difficult. The poor statistics available for these properties means that their determination requires either very long equilibrium simulations, or the use of nonequilibrium molecular dynamics. Experimentally, the shear viscosity of water has been well studied. In contrast, although many models have been proposed for the study of liquid water, which reproduce a number of properties of liquid water, the viscosity of these models has not (to our knowledge) been determined. This is due, in part, to the poor statistics inherent to these methods, and to
the confusion concerning the most efficient method for the determination of shear viscosity from an equilibrium simulation [ 2-5 1. Here we report the determination of the zero-frequency shear viscosity for the SPC [ 61 and SPC/E [ 71 models of water at two different temperatures, 277 and 300 K, using long equilibrium molecular dynamics simulations. The effect of including long-range interactions using the reaction field technique was also investigated.
2. Method A total of eight molecular dynamics simulations were performed. The SPC [6] and SPC/E [ 71 models of water were simulated at 277 and 300 K using a reaction field potential [ 81 with a permittivity of either eRr= 1 (no reaction field) and eRF= cc. All simulations involved a system of 5 12 waters in a cubic periodic box. The systems were initially equilibrated in the NPT ensemble for 50 ps. The average box lengths from these simulations were then used as the box dimensions for a series of NVT simulations. A further 50 ps equilibration was performed in the NVT ensemble. Each system was then simulated for a further 1000 ps. All simulations were performed using the GROMOS program [ 9 1, with a cutoff distance of Rc = 0.9nm. Temperature regulation was achieved with the weak coupling method [ lo] and a coupling constant of 0.1 ps. Bond lengths and angles were constrained using SHARE [ 111 and
0009-261419316 06.00 0 1993 Elsevier Science Publishers B.V. All rights reserved.
315
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
a relative tolerance of 10m5.The time step was 2 fs. The shear viscosity is related to the rate of momentum transport through the system. The shear viscosity can be obtained from simulation by either Green-Kubo- or Einstein-type expressions. The Green-Kubo approach relates the shear viscosity q to the time autocorrelation function of the off-diagonal elements of the pressure tensor [ 1,12 1,
3 December 1993
t
AG(t)=G(t)-G(O)=
j”i;(t’) dt’,
(5)
0
where the time derivative of G(t) is given by .
G(t)=
’
7
=Pa/3(0 , (1)
where Y is the volume of the system, kB is the Boltzmann constant, T is the absolute temperature, t is time, and Pap is the a/3 component of the pressure tensor (cr,+x, y, z; cr#B) given by
PC@(t) = _: x
5 ( i=l
+ mi
“d ‘= f.“, 2k,Tdt
5 icj
Gja(wj/.d~)
3
(2)
>
(
[G(t)-G(0)12) ,
(3)
where G( t ) is given by
ijglria(t)Pi,dt)
and ri( t) and pi( t) are the position and momentum of the centre of mass of molecule i at time t, respectively. The angular brackets in eqs. ( 1) and (3) denote an average over time origins. Unfortunately, G(t) is not continuous under periodic boundary conditions. However, one may write [ 161 316
and& is the total force on atom i. Then one has an Einstein-type relationship which is independent under periodic boundary conditions. This is the approach we have adopted here. During the simulations, each element of the pressure tensor was saved at every step for analysis. We have used these elements of the pressure tensor to calculate the viscosity using eqs. (3), (5) and (6).
Pia(f>Pi#(t>
pr( t) is the centre of mass momentum of molecule i (mass mi) at time t, rij= ri - rj is the intermolecular distance vector, and& is the force on molecule i by molecule j, and the summation is over all N molecules. However, the above correlation function often contains significant long time tails [ 13 1, which have to be integrated correctly to obtain a reliable value of ?/. To avoid the problem of long time tails many transport properties are determined using the equivalent Einstein expression. Helfand derived the following, generally accepted [ 35 1, Einstein relationship for the shear viscosity [ 14,15 ] :
G(t)=
(6)
3. Results The time development of AG(t)’ from one of the simulations is shown in fig. 1. The other simulations produced data of a similar quality. After an initial nonlinear regime of about 5 ps, the curves become linear in time. The shear viscosity was calculated 2000 ,
I
I
J
1500 ?-.
a
“E
s
P 1000 a n “E
s? ”
500
0
0
2
6
4 Time
8
10
(ps)
Fig. 1. The mean-square displacement of G(t) as a function of time obtained by integrating the off-diagonal elements of the pressure tensor. The curves are for the SPC/E model simulated at 300 K without reaction field. The thin lines correspond to the three independent off-diagonal elements of the pressure tensor, and the thick line is the resulting average.
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
from the slope of these curves obtained between 5 10 ps. Errors were determined from the variation between the results for the different off-diagonal components. The final values are presented in table 1. The results for the SPC models are not so good. The viscosity is too low at both temperatures and there is too small a variation of viscosity with temperature. The SPC/E model, however, performs very well with excellent agreement at 300 K and acceptable agreement at 277 K, although the viscosity is slightly too low. The above results reflect nicely the differences between the SPC and SPC/E water models. The SPC/ E model has slightly larger partial atomic charges than the SPC model (qo= -0.8476 versus -0.8200 e). Hence, the SPC/E model has a higher density, a lower diffusion constant and a larger viscosity than the SPC model. For all three properties the SPC/E model gives better agreement with experiment. It is interesting to note that using a reaction field treatment of the long-range interactions did not appear to have any significant effect on the viscosity. However, the reaction field did affect the self-diffusion constant of water as observed in the simulations (see table 1). It is known that there is an in-
and
timate relationship [ 19,201 between viscosity and self-diffusion (Da 1 /q ) , and therefore one might expect a significant variation in viscosity with eRF.Table 1 shows the product Dr7u/k,J, where d is a molecular diameter. This product should be constant if Stokes law is obeyed [ 191. Experimentally, it is constant over the temperature range covered here. However, the simulation values vary by as much as 30% between different models and reaction fields, although this is probably due to the sizeable errors estimated for the calculations. In fig. 2 we display the pressure autocorrelation functions calculated for the SPC/E model at 277 and 300 K without a reaction field treatment. Both the diagonal and off-diagonal components of the pressure tensor are shown in fig. 2a. The diagonal elements display a large initial loss in correlation, followed by a very slow decay over many picoseconds. The off-diagonal elements display an even faster initial decay and a smaller long-time correlation. Both types of correlation function contain some periodicity; a periodicity common to both functions of about 0.10 ps, and a smaller period for the off-diagonal elements of about 0.05 ps. In fig. 2b the difference in the correlation functions with tempera-
Table 1 Molecular dynamics results
T(K)
fIlF
Box length L
Density
Diffusion m*/s)
SPC 277 300
1 co 1 W
SPC/E 277
1 W
300
exp. 0) 277 300
1 co
Viscosity
Pressure
(cp) *)
(W/mol/nm3)
Dqu b)'kBT
P
2.49015 2.50756 2.49975 2.52286
0.991 0.971 0.980 0.953
2.80 3.59 3.94 5.28
0.79( 16) 0.80(4) 0.54(9) 0.58(9)
-3.27 0.85 1.73 0.11
0.223 0.292 0.199 0.289
2.47209 2.49107 2.48133 2.50257
1.013 0.990 1.002 0.976
1.62 1.75 2.68 3.18
1.39(6) 1.40(21) 0.82(9) 0.91(7)
-1.65 -0.28 -0.28 -2.26
0.226 0.247 0.204 0.271
1.000 0.996
1.1 *I 2.2 *)
1.57 0.85
0.174 0.174
a) Standard deviations are given in parentheses. b, The molecular diameter u was obtained from 2 x(IY/~)~= L 3/5 12. Cl Ref. 1171. *‘Ref. 1181.
317
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
3 December 1993
found to be negligible. The results suggest that the SPC model should not be used for low-temperature studies, or for simulations where hydrodynamical effects are critical.
2000
1500
1000
Acknowledgement
500
0
-500 t 0.0
0.2
0.4
0.6
0.6
1.0
1000
We thank Ilario G. Tironi for valuable discussions and gratefully acknowledge the generous financial support of Zeneca Ltd which made this study possible.
600
600
References 400
[ 1] R. Zwanzig, Ann. Rev. Phys. Chem. 16 (1965) 67. [ 21 D.A. McQuarrie, Statistical mechanics (Harper and Row,
200
0
.200' 0.0
a 0.2
s 0.4
' 0.6
s 0.6
' 1.0
Time(ps)
Fig. 2. Pressure time autocorrelation functions for the SPC/E model simulated without reaction field. In (a) the upper curves correspond to the diagonal elements (xx, yy and zz), and the lower curves correspond to the off-diagonal elements (xy, xz and yz) at 300 K. In (b) the thick lines correspond to 300 K and the thin lines to 277 K. The top two curves correspond to the xx component of the pressure tensor and the bottom two curves to the xy component.
ture is displayed. Differences in the correlation functions appear at long times and are greater for the diagonal elements. Integration of the off-diagonal elements to 10 ps gave viscosities which were essentially within the error of those calculated using the Einstein-type relationship. Integration of the diagonal elements (to give the bulk viscosity) was not possible.
4. Conclusions The zero-frequency shear viscosity of SPC and SPC/E water has been determined using long equilibrium molecular dynamics simulations. The SPC/ E model reproduces the experimental viscosity at both 277 and 300 IS. The SPC model does not. The effect of long-range interaction on the viscosity was 318
New York, 1976). [3]A.A.ChialvoandP.G.Debenedetti,Phys.Rev.A43 (1991) 4289. [4] A.A. Chialvo, P.T. Cummings and D.J. Evans, Phys. Rev. E 47 (1993) 1702. [5]M.P.AllenandA.J. Masters,Mol.Phys. 79 (1993) 435. [6] H.J.C. Beret&en, J.P.M. Postma, W.F.vanGunsterenand J. Hermans, in: Intermolecular forces. Interaction models for water in relation to protein hydration, ed. B. Pullman (Reidel, Dordrecht, 1981) pp. 331-342. [7] H.J.C. Berendsen, J.R. GrigeraandT.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [8] J.A. BarkerandR.0. Watts, Mol. Phys. 26 (1973) 789. [9] W.F. van Gunsteren and H.J.C. Berendsen, Groningen molecular simulation (GRO-MOS) Library Manual, Biomos, Groningen (1987). [ lo] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola and J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [ 111 J.-P. Ryckaert, G. Ciccotti and H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. 1121 F.C. Andrews, J. Chem. Phys. 47 (1967) 3161. [ 131 B.J. Alder, in: Molecular dynamics simulation of statistical mechanical systems. Molecular dynamics simulations, eds. G. Ciccotti and W.G. Hoover (North-Holland Amsterdam, 1986) pp. 66-80. [14]E.Helfand,Phys.Rev. 119 (1960) 1. [ 151 D. Gass, J. Chem. Phys. 51 (1969) 4560. [ 16 ] B.J. Alder, D.M. Gass and T.E. Wainwright, J. Chem. Phys. 53 (1970) 3813. [ 171 R.C. Weast, Handbook of chemistry and physics (CRC Press, Boca Raton, 1975). [ 18 ] K.T. Gillen, DC. Douglas and M.J.R. Ho&, J. Chem. Phys. 57 (1972) 5117. [ 191 R. Zwanzig, J. Chem. Phys. 79 (1983) 4507. [20] R. Zwanzig and AK. Harrison, J. Chem. Phys. 83 (1985) 5861.