Advanced
Temperature Dependence on Structure and Self-Diffusion of Water: A Molecular Dynamics Simulation Study using SPC/E Model
Temperature Dependence on Structure and Self-Diffusion of Water: A Molecular Dynamics Simulation Study using SPC/E Model
Bulletin of the Korean Chemical Society. 2013. Dec, 34(12): 3800-3804
Copyright © 2013, Korea Chemical Society
  • Received : August 26, 2013
  • Accepted : October 02, 2013
  • Published : December 20, 2013
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Song Hi Lee

Abstract
In this study, molecular dynamics simulations of SPC/E (extended simple point charge) model have been carried out in the canonical NVT ensemble over the range of temperatures 300 to 550 K with and without Ewald summation. The quaternion method was used for the rotational motion of the rigid water molecule. Radial distribution functions g OO (r), g OH (r), and g HH (r) and self-diffusion coefficients D for SPC/E water were determined at 300-550 K and compared to experimental data. The temperature dependence on the structural and diffusion properties of SPC/E water was discussed.
Keywords
Introduction
Liquid water is the most important solvent in nature that has many special and unusual physicochemical properties. Most of these water properties are due to the ability of water molecules to form hydrogen bonds to other water molecules in three-dimensional directions. In the past decades, many classical force fields for molecular simulations on water have been developed. 1 - 13 A 2002 review indicates that there are 46 water models, 14 which were classified as rigid, flexible and polarizable models. 15 The most popular water models - the TIP3P (transferable intermolecular potential 3P) (original 1 and modified 2 ), SPC (simple point charge) (original 3 and refined 4 ), and SPC/E (extended simple point charge) 5 can be described as effective rigid pair potentials composed of Lennard-Jones and Coulombic terms.
Recently, many molecular dynamics (MD) studies for selfdiffusion coefficient, 16 - 23 viscosity, 23 - 31 and thermal conductivity 31 - 34 using various water models have been reported. Our interest in this study concentrated on the radial distribution functions, g OO (r), g OH (r), and g HH (r) and self-diffusion coefficient D of the bulk water. The experimental values for self-diffusion coefficients of pure water have been measured to be between 2.26 and 2.29 (× 10 −9 m 2 /s) 35 - 37 at 298.15 K. Self-diffusion coefficients have been reported using MD simulations for the original TIP3P water model between 5.2 and 7.0 (× 10 −9 m 2 /s), 16 for the modified TIP3P water model between 2.3 and 5.2 (× 10 −9 m 2 /s), 17 - 21 for the original SPC water model between 3.6 and 5.2 (× 10 −9 m 2 /s), 16 for the refined SPC water model between 4.2 and 4.4 (× 10 −9 m 2 / s), 22 and for SPC/E water model between 2.2 and 4.4 (× 10 −9 m 2 /s) 16 at 298 K. However, the temperature dependence of self-diffusion coefficient at high temperatures is hardly found in the literature except Ref.38 over the range of temperatures 273 to 373 K, even though MD simulation studies for selfdiffusion coefficient in supercritical water 39 40 have been reported.
In this study, we utilize the Einstein and Green-Kubo relations for the calculation of self-diffusion coefficients of SPC/E water using MD simulations over the range of temperatures 300 to 550 K. The primary goal of this study is to compare self-diffusion coefficients of water with the experimental measures at high temperatures and to examine the temperature dependence of the radial distribution functions and self-diffusion coefficients of SPC/E water. We describe the molecular models and the technical details of MD simulation in section Ⅱ, our results in section Ⅲ, and the concluding remarks in section Ⅳ.
Molecular Models and Molecular Dynamics Simulation
Water is simulated using the extended simple point charge (SPC/E) model 5 in which the charges on H are at 1.000 Å from the Lennard-Jones center at O, the negative charge is at the O site, and the HOH angle is 109.47°. The pair potential between water and ion has the form
PPT Slide
Lager Image
where σ = 3.169 Å and ε = 0.6502 kJ/mol are Lennard-Jones (LJ) parameters between oxygen atoms on different water molecules, q i is the charge at site i in water ( q O = −0.8476 e and q H = 0.4238 e ), and r ij is the distance between charge sites i and j in different water molecules.
The preliminary canonical ensemble (NVT fixed) MD simulations of N = 1024 water molecules with and without Ewald summation over the range of temperatures 300 to 550 K were started for equilibration in the cubic box of length L determined from water densities at given temperatures (see Table 1 ). Ewald summations were used in our simulations with the parameter for κ = 2.0 Å −1 and the real-space cut distance r cut and K max chosen as 10.0 Å and 7, respectively. Nose-Hoover thermostat 41 42 was used to keep the temperature constant (the Nose-Hoover thermostat relaxation constant is given as Q = 10 f k B with f as the number of degrees of freedom and k B as the Boltzmann constant). The usual periodic boundary condition was applied in the x-, y-, and z-directions, and the minimum image convention for pair potential were applied. The equations of translational motion were solved using a velocity Verlet algorithm 43 with a time step of 10 −15 second (1 fs) and a quaternion formulation 44 45 was employed to solve the equations of rotational motion about the center of mass of rigid SPC/E water molecules. The configurations of water molecules were stored every 10 time steps for further analysis. The systems were fully equilibrated for 500,000 time steps and the equilibrium properties are averaged over 10 blocks of 100,000 time steps (0.1 ns).
Densities (ρm, kg/m3) and lengths (L, Å) of cubic simulation box of each system at given temperatures T(K), and comparison of self-diffusion coefficients (10−9m2/s) calculated from MD simulations with and without Ewald summation to experimental measures. Uncertainties (standard deviation) in the last reported digit(s) are given in the parenthesis
PPT Slide
Lager Image
aRefs.37 and 48. bRefs. 37, 48, and 49.
There are two routes to determine self-diffusion constants of water from MD simulations; the Einstein relation from the mean square displacement (MSD), 46
PPT Slide
Lager Image
and the Green-Kubo relation from the velocity autocorrelation (VAC) function. 46
PPT Slide
Lager Image
Results and Discussion
The radial distribution functions, g OO (r), g OH (r), and g HH (r) computed from our MD simulations for the SPC/E water at 300, 400, and 500 K are compared with neutron and x-ray diffraction data 47 at 298 and 423 K in Figures 1 , 2 , and 3 , respectively. Heights and positions of the peaks and the first minima at 300-550 K are given for g OO and g OH in Table 2 . The overall trend is that the peak heights are lowered and the valley heights are raised for all the three radial distribution functions with increasing temperature. The positions of the first and second maxima and the first minima for g OO (r) and g OH (r) are shifted to longer distances with increasing temperature except those of the first minima of g OH (r) which are never changed with temperature, as seen in Table 2 .
Positions and Magnitudes at Maxima and Minima of gOOand gOHradial distribution functions at 300-550 K using SPC/E model
PPT Slide
Lager Image
Positions and Magnitudes at Maxima and Minima of gOO and gOH radial distribution functions at 300-550 K using SPC/E model
When compared with neutron and x-ray diffraction data, 47 the MD result of g OO (r) at 300 K in Figure 1 shows an overall good agreement with the experiment at 298 K except the first peak is slightly higher, while the agreement between the MD result at 400 K and the experiment at 423 K is poor with shorter position and lower height of the first peak in the MD g OO (r). For g OH (r) in Figure 2 , the MD result at 300 K has a too high first peak but the first minimum and the second peak are acceptable with the experiment at 298 K. The agreement between the MD result at 400 K and the experiment at 423 K is also poor with a shorter position and a higher height of the first peak and a lower height of the first minimum in the MD g OH (r). The agreement between MD result of g HH (r) at 300 K and the experiment at 298 K is in a good accordance as seen in Figure 3 , while that between the MD at 400 K and the experiment at 423 K is also poor, and especially the deviation for the first minimum height is the worst. The overall disagreement between the MD result at 400 K and the experiment at 423 K for all the three radial distribution functions may not be attributed to the temperature difference.
PPT Slide
Lager Image
MD results of gOO(r) at 300-500 K and the experiment values at 298 and 423 K.
PPT Slide
Lager Image
MD results of gOH(r) at 300-500 K and the experiment values at 298 and 423 K.
Self-diffusion coefficients evaluated from the slopes of mean square displacements (MSD) of SPC/E water using Eq. (2) and the time-integrations of velocity auto-correlation (VAC) functions using Eq. (3) at 300-550 K with and without Ewald summation are shown in Table 1 . The calculated MSDs ( Fig 4 .) with Ewald summation show a perfect linear behavior after 0.5 ps and D was obtained from the slopes of MSDs between 2 and 10 ps, and also the tails of VACs ( Fig 5 .) with Ewald summation decay to zero very quickly within 1 ps and D was obtained from the averages of the timeintegrations of VACs from 0 to 2 ps and from 0 continually to10 ps. The behaviors of MSDs and VACs (not shown) without Ewald summation are very similar to those with Ewald summation. Four values of D at 300-550 K from MSDs and VACs with and without Ewald summation are almost the same except that D obtained with Ewlad summation are slightly larger than those without it as seen in Table 1 .
PPT Slide
Lager Image
MD results of gHH(r) at 300-500 K and the experiment values at 298 and 423 K.
PPT Slide
Lager Image
Mean square displacements (MSD) of SPC/E water at 300-550 K.
PPT Slide
Lager Image
Velocity auto-correlation (VAC) functions of SPC/E water at 300-550 K.
Two sets of experimental measures 37 48 49 for D of water are also listed in Table 1 . The fitted function for D of water as a function of temperature T(K) was of the form 37 :
PPT Slide
Lager Image
where a 0 = 6.11903, a 1 = −1.195593, a 2 = −0.05229, and a 3 = −0.01841 for D(Exp1) over 278-498 K and a 0 = 13.2172, a 1 = −9.08602, a 2 = 2.80883, and a 3 = −0.35713 for D(Exp2) over 242-498 K.
PPT Slide
Lager Image
Comparison of diffusion coefficients (D) with experimental data as a function of inverse temperature.
In Figure 6 , we compared D obtained from our MD simulations with the experimental data as a function of inverse temperature. The overall MD results underestimates D except at 300 K and the deviation increases with increasing temperature, 5% (350 K), 17% (400 K), 24% (450 K), and 25% (500 K), when compared D(MSD) with Ewald summation to D(Exp1) in Table 1 . This behavior of D with temperature is in a good accord with the previous study over the range of temperatures 273 to 373 K, 38 and especially the cross point between the MD result for D and the experiment D coincides exactly at T = 330 K 38 and at 1000/T = 3 in Figure 3 . However the reason for the deviation for D at high temperatures is not fully understood. Perhaps the potential parameters for the SPC/E model might be developed at room temperature.
The discrepancy between self-diffusion coefficients of SPC/ E water and the experiment measures at high temperatures requires a more refined water model for self-diffusion coefficients. A recent study 23 using a rigid non-polarizable water model, TIP4P/2005, 50 for diffusion coefficient D and shear viscosity ƞ of rigid water models has reported an excellent agreement with the experimental results at 300 K for D and ƞ. MD simulations for D and ƞ using TIP4P/2005 at high temperatures are currently under investigation.
The temperature dependence of the MD result for D and the experiment D over the range of temperatures 300 to 500 K are suitably described by an Arrhenius plot as shown in Figure 6 :
PPT Slide
Lager Image
where D 0 is the pre-exponential factor, RT has the usual meaning, and E D is the activation energy of water diffusion. The value of the activation energy is a direct measure of the temperature dependence of self-diffusion coefficient. The calculated activation energies are 3.1 and 3.7 kcal/mol for D(MSD) with Ewald summation and for D(Exp1), respectively. The activation energies reported for the experiment D of water over the range of temperatures 273 to 323 K 36 37 51 and for over the range of temperatures 243 to 298 K 49 are 4.3 and 1.6 kcal/mol, respectively.
Conclusion
We have carried out molecular dynamics of SPC/E model in the canonical NVT ensemble over the range of temperatures 300 to 550 K with and without Ewald summation. The overall trend of the calculated radial distribution functions, g OO (r), g OH (r), and g HH (r) is that the peak heights are lowered and the valley heights are raised with increasing temperature. When compared with experimental data, the MD results for all the three radial distribution functions at 300 K show an overall good agreement with the experiment at 298 and the overall disagreement between the MD result at 400 K and the experiment at 423 K may not be attributed to the temperature difference.
Self-diffusion coefficients evaluated from the slopes of mean square displacements (MSD) and the time-integrations of velocity auto-correlation (VAC) functions at 300-550 K with Ewald summation are compared with the experimental data. The overall MD results underestimate D except at 300 K and the deviation increases with increasing temperature. The temperature dependence of the MD result for D and the experiment D over the range of temperatures 300 to 500 K are suitably described by an Arrhenius plot. The calculated activation energies are 3.1 and 3.7 kcal/mol for D(MSD) with Ewald summation and for D(Exp1), respectively.
Acknowledgements
This research was supported by Kyungsung University Research Grants in 2013.
References
Jorgensen W. L. , Chandrasekhar J. , Madura J. D. , Impey R. W. , Klein M. L. 1983 J. Chem. Phys. 79 926 -    DOI : 10.1063/1.445869
Neria E. , Fischer S. , Karplus M. 1996 J. Chem. Phys. 105 1902 -    DOI : 10.1063/1.472061
Berendsen H. J. C. , Postma J. P. M. , van Gunsteren W. F. , Hermans J. 1981 In Intermolecular Forces; Pullman, B., Ed. Reidel Dordrecht
Berweger C. D. , van Gunsteren W. F. , Müller-Plape F. 1995 Chem. Phys. Lett. 232 429 -    DOI : 10.1016/0009-2614(94)01391-8
Berendsen H. J. C. , Grigera J. R. , Straatsma T. P. 1987 J. Phys. Chem. 91 6269 -    DOI : 10.1021/j100308a038
1976 J. Chem. Phys. 64 1351 -    DOI : 10.1063/1.432402
Watanabe K. , Klein M. L. 1989 Chem. Phys. 131 157 -    DOI : 10.1016/0301-0104(89)80166-1
Liu Y. , Ichiye T. 1996 J. Phys. Chem. 100 2723 -    DOI : 10.1021/jp952324t
Buch V. , Sandler P. , Sadlej J. 1998 J. Phys. Chem. B 102 8641 -
Levitt M. , Hirshberg M. , Sharon R. , Laidig K. E. , Daggett V. 1997 J. Phys. Chem. B 101 5051 -    DOI : 10.1021/jp964020s
Jorgensen W. L. , Jenson C. 1998 J. Comput. Chem. 19 1179 -    DOI : 10.1002/(SICI)1096-987X(19980730)19:10<1179::AID-JCC6>3.0.CO;2-J
Chialvo A. A. , Cummings P. T. 1996 J. Chem. Phys. 105 8274 -    DOI : 10.1063/1.472718
Dang L. X. 1998 J. Phys. Chem. B 102 620 -
Guillot B. 2002 J. Mol. Liq. 101 219 -    DOI : 10.1016/S0167-7322(02)00094-6
Caleman C. 2007 J. Chem. Phys. 22 709 -
van der Spoel D. , van Maaren P. J. , Berendsen H. J. C. 1998 J. Chem. Phys. 108 10220 -    DOI : 10.1063/1.476482
Feller S. E. , Pastor R. W. , Rojnuckarin A. , Bogusz S. , Brooks B. R. 1996 J. Phys. Chem. 100 17011 -    DOI : 10.1021/jp9614658
Tasaki K. , McDonald S. , Brady J. W. 1993 J. Comput. Chem. 14 278 -    DOI : 10.1002/jcc.540140304
Liu Q. , Schmidt R. K. , Teo B. , Karplus P. A. , Brady J. W. 1997 J. Am. Chem. Soc. 119 7851 -    DOI : 10.1021/ja970798v
Smip P. E. , Blatt H. D. , Pettitt B. M. 1997 J. Phys. Chem. B 101 3886 -    DOI : 10.1021/jp9637643
Makarov V. A. , Feig M. , Andrews B. K. , Pettitt B. M. 1998 Biophys. J. 75 150 -    DOI : 10.1016/S0006-3495(98)77502-2
Mark P. , Nilsson L. 2001 J. Phys. Chem. A 105 9954 -    DOI : 10.1021/jp003020w
Tazi S. , Bopan A. , Salanne M. , Marry V. , Turq P. , Rotenberg B. 2012 J. Phys.: Condens. Matter 24 284117 -    DOI : 10.1088/0953-8984/24/28/284117
Balasubramanian S. , Mundy C. J. , Klein M. L. 1996 J. Chem. Phys. 105 11190 -    DOI : 10.1063/1.472918
Zhang Y.-G. 2001 Mol. Phys. 99 283 -    DOI : 10.1080/00268970010011762
Hess B. 2002 J. Chem. Phys. 116 209 -    DOI : 10.1063/1.1421362
Wu Y. J. , Tepper H. L. , Vop G. A. 2006 J. Chem. Phys. 124 024503 -    DOI : 10.1063/1.2136877
Chen T. , Smit B. , Bell A. T. 2009 J. Chem. Phys. 131 246101 -    DOI : 10.1063/1.3274802
Wensink E. , Hoffmann A. , van Maaren P. , van der Spoel D. 2003 J. Chem. Phys. 119 7308 -    DOI : 10.1063/1.1607918
González M. A. , Abascal J. L. F. 2010 J. Chem. Phys. 132 096101 -    DOI : 10.1063/1.3330544
Zhang Y. 2012 Chem. Phys. Lett. 542 37 -    DOI : 10.1016/j.cplett.2012.05.044
Tani A. 1995 Phys. Rev. E 51 1091 -    DOI : 10.1103/PhysRevE.51.1091
Isachenko V. P. , Osipova V. A. , Sukomel A. S. 1980 Heat Transfer Mir Moscow
Kataoka Y. 1989 Bull. Chem. Soc. Jpn. 62 1421 -    DOI : 10.1246/bcsj.62.1421
Atkins P. , Paula J. D. 2002 Physical Chemistry 7th ed Freeman New York
Holz M. , Heil S. R. , Sacco A. 2002 Phys. Chem. Chem. Phys. 2 4740 -
Easteal A. J. , Price W. E. , Woolf L. A. 1989 J. Chem. Soc., Faraday Tans. 1 85 1091 -    DOI : 10.1039/f19898501091
Peris C. S. , Coorey R. V. , Weerasinghe S. 2010 Proc. Tech. Sess. 26 75 -
Yoshida K. , Matubayasi N. , Uosaki1 Y. , Nakahara M. 2010 J. Phys.: Condens. Matter 215 012093 -
Lee S. H. 2013 Bull. Kor. Chem. Soc. 36
Nosé S. 1984 J. Chem. Phys. 81 511 -    DOI : 10.1063/1.447334
Hoover W. G. 1985 Phys. Rev. A 31 1695 -    DOI : 10.1103/PhysRevA.31.1695
Swope W. C. , Andersen H. C. , Beren P. H. , Wilson K. R. 1982 J. Chem. Phys. 76 637 -    DOI : 10.1063/1.442716
Evans D. J. 1977 Mol. Phys. 34 317 -    DOI : 10.1080/00268977700101751
Evans D. J. , Murad S. 1977 Mol. Phys. 34 327 -    DOI : 10.1080/00268977700101761
McQuarrie D. A. 1976 Statistical Mechanics Harper and Row New York
Soper A. K. 2000 Chem. Phys. 258 121 -    DOI : 10.1016/S0301-0104(00)00179-8
Kell G. S. 1977 J. Phys. Chem. Ref. Data 6 1109 -    DOI : 10.1063/1.555561
Gillen K. T. , Douglass D. C. , Hoch M. J. R. 1972 J. Chem. Phys. 57 5117 -    DOI : 10.1063/1.1678198
Abascal J. L. F. , Vega C. 2005 J. Chem. Phys. 123 234505 -    DOI : 10.1063/1.2121687
Lee S. H. , Rasaiah J. C. 2011 J. Chem. Phys. 135 124505 -    DOI : 10.1063/1.3632990