Advanced
Molecular Dynamics Simulation Study for Shear Viscosity of Water at High Temperatures using SPC/E Water Model
Molecular Dynamics Simulation Study for Shear Viscosity of Water at High Temperatures using SPC/E Water Model
Bulletin of the Korean Chemical Society. 2014. Feb, 35(2): 644-646
Copyright © 2014, Korea Chemical Society
  • Received : October 13, 2013
  • Accepted : November 17, 2013
  • Published : February 20, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Song Hi Lee

Abstract
Keywords
Molecular Dynamics Simulation
Canonical ensemble MD simulations of N = 1024 water molecules with Ewald summation over the range of temper-atures 300 to 550 K were carried out 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. Nosé-Hoover thermostat 30 31 was used to keep the temperature constant (the Nosé-Hoover thermostat relaxation constant is given as Q = 10 f kB with f as the number of degrees of freedom and kB as the Boltzmann constant). The usual periodic boundary condition was applied in the x-, y-, and z-direction, and the minimum image convention for pair potential were applied. The equations of translational motion were solved using a velocity Verlet algorithm 32 with a time step of l0 -15 second (1 fs) and a quaternion formu-lation 33 34 was employed to solve the equations of rotational motion about the center of mass of rigid SPC/E water mole-cules. 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 viscosity (cP) calculated from MD simulations with Ewald summation to experimental measures. Uncertainties (standard deviation) in the last reported digit(s) are given in the parenthesis
PPT Slide
Lager Image
aRef.22
The original GK formula for viscosity ƞ is given by
PPT Slide
Lager Image
where
PPT Slide
Lager Image
is the pressure of particle i on to the wall with αβ = xy, xz, yx, yz, zx, and zy.
Results and Discussion
We plot the stress auto-correlation (SAC) functions, the integrand of Eq. (1), of SPC/E water at very short times for the temperatures of 300, 400, and 500 K in Figure 1 . The SAC functions have lots of oscillations with a fast decay within 0.06 ps and an oscillatory peak at 0.13 ps. These functions are somewhat similar to the ones reported by Guo et al . 27 and Tazi et al . 29 which show much more smooth curve without lots of oscillations. The positions of the first minimum (at 0.06 ps) and the first peak (at 0.13 ps) coincide exactly. At 300 K, the SAC function does not decay to zero at 0.8 ps which reflects a large value of shear viscosity. As the temperature increases, the SAC functions lowered with-out changing the magnitude and the position of the first minimum.
PPT Slide
Lager Image
Stress auto-correlation (SAC, (kJ/mol·nm3)2) functions of SPC/E water at 300-500 K. The inset shows the long-time behavior of SAC functions.
Running integrals for ƞ(t) of SPC/E water at 300-550 K are plotted as a function of time in Figure 2 . All the running integrals for viscosity clearly show plateaus which signify that the corresponding SAC functions have decayed to zero and are fluctuating along the horizontal time axis. In general, ƞ(t) for longer correlation will have larger statistical un-certainty because less data are available for its calculation. Therefore, the beginning of a plateau, which corresponds to the time when the SAC function reaches zero (not the two zero values before 0.2 ps), represents the start of the desired value of viscosity with the smallest uncertainty. As shown in the inset of Figure 1 , all the SAC functions reach zero at about 4 ps and we report the shear viscosities at 300-550 K in Table 1 by averaging the running integrals for viscosity in Figure 2 for 4-8 ps.
PPT Slide
Lager Image
Running integrals for ƞ(t) of SPC/E water at 300-550 K.
The shear viscosities obtained by MD simulations at room temperature are discussed above. Our result for the shear viscosity using SPC/E potential at 300 K, 0.722 cP, differs from the MD results using other potentials – TIP3P, TIP4P, TIP5p, and TIP4P/2005, but agrees well with the MD results using SPC/E. It is very interesting to compare the shear viscosities evaluated using various water potentials, all the values having very similar values for each water potential.
We have compared our MD results for viscosity ƞ with experimental data as a function of inverse temperature in Figure 3 . The shear viscosity underestimates the experimental measures at low temperatures and overestimates at high temperatures; (ƞ MD − ƞ exp )/ƞ exp = −15% (300 K), −8% (350 K), 15% (400 K), 16% (450 K), 20% (500 K), and 5% (550 K). The temperature dependence of the MD result for ƞ MD and the experiment ƞ Exp over the range of temperature 300 to 550 K are suitably described by an Arrhenius plot as shown in Figure 3 : ƞ = ƞ o exp(E ƞ /RT) where ƞ o is the pre-exponential factor and E ƞ is the activation energy of water viscosity. The value of the activation energy is a direct mea-sure of the temperature dependence of shear viscosity. The calculated activation energies are 2.5 and 2.9 kcal/mol for ƞ MD with Ewald summation and for ƞ Exp , respectively. The corresponding diffusion activation energies [D = D o exp(−E D /RT)] are 3.1 and 3.7 kcal/mol for D MD and for D Exp , respectively. 35
According to the Stokes-Einstein (SE) relation about a Brownian (B) particle immersed in a solvent (s), D B = k B T/ Cπƞ s r B , where C is the hydrodynamic boundary condition, and ƞ s and r B are the viscosity of solvent and the radius of the B particle. Applying the SE relation to bulk water, ƞD ∝ k B T. We plot ƞD vs T in the inset of Figure 3 . The behavior of ƞD for the experimental data shows a linear increment with T, but that for the MD result does not, which reflects the requirement of a more refined water model for self-diffusion coefficients and shear viscosities at high temperatures. MD simulations for D and ƞ using TIP4P/2005 16 at high temper-atures are currently under investigation.
PPT Slide
Lager Image
Comparison of viscosity ƞ with experimental data as a function of inverse temperature. The error bar indicates standard deviation. The inset shows ƞD (10-9 cP·m2/s) as a function of T.
In contrast for the Green-Kubo (GK) formula to predict single-particle properties such as self-diffusion coefficient through the velocity auto-correlation (VAC) function as accurately as N (number of particle) times, the GK formula seems to have difficulty in predicting collective properties such as shear viscosity and thermal conductivity through the stress (SAC) and heat flux auto-correlation (HFAC) functions since the collective properties are the system properties which have N times less accurate statistics when compared with the single-particle properties. From this point of view, while the large deviation of our previous result for self-diffusion coefficients at 300-550 K 35 using SPC/E from the experimental data by -25~16% may come from the SPC/E model not from the GK formula, the current result for shear viscosities at the same temperatures using the same water model deviates from the experimental data by -25~20% which is rather a better agreement than those for self-diffusion coefficient since the stress is a collective property. Alternative way to overcome the poor statistical accuracy in evaluating shear viscosity and thermal conductivity through the stress (SAC) and heat flux auto-correlation (HFAC) functions is to consider the stress and heat flux as single-particle properties not as collective properties. Then the statistical accuracy is improved by N (number of particle) times large. The modified GK formula for the better statistical accuracy was applied to the Lennard-Jones systems. 36 37
In summary, the shear viscosities of SPC/E water model over the range of temperatures 300 to 550 K were evaluated using the Green-Kubo formula by carrying out molecular dynamics (MD) simulations in the canonical NVT ensemble with Ewald summation. The MD result for the shear visco-sity underestimates the experimental measures at low temper-atures and overestimates at high temperatures. The temper-ature dependence of the MD result for ƞ MD and the experi-ment ƞ Exp over the range of temperature 300 to 550 K are suitably described by an Arrhenius plot. The calculated activation energies are 2.5 and 2.9 kcal/mol for ƞ MD with Ewald summation and for ƞ Exp , respectively.
References
Jorgensen W. L. , Chandrasekhar J. , Madura J. D. , Impey R. W. , Klein M. L. 1983 J. Chem. Phys. 79 926 -
Neria E. , Fischer S. , Karplus M. 1996 J. Chem. Phys. 105 1902 -
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-Plathe F. 1995 Chem. Phys. Lett. 232 429 -
Berendsen H. J. C. , Grigera J. R. , Straatsma T. P. 1987 J. Phys. Chem. 91 6269 -
Matsuoka O. , Clementi E. , Yoshimine M. 1976 J. Chem. Phys. 64 1351 -
Watanabe K. , Klein M. L. 1989 Chem. Phys. 131 157 -
Liu Y. , Ichiye T. 1996 J. Phys. Chem. 100 2723 -
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 -
Jorgensen W. L. , Jenson C. 1998 J. Comput. Chem. 19 1179 -
Chialvo A. A. , Cummings P. T. 1996 J. Chem. Phys. 105 8274 -
Dang L. X. 1998 J. Phys. Chem. B 102 620 -
Guillot B. 2002 J. Mol. Liq. 101 219 -
Caleman C. 2007 J. Chem. Phys. 22 709 -
Abascal J. L. F. , Vega C. 2005 J. Chem. Phys. 123 234505 -
Vega C. , Abascal J. L. F. , Nezbeda I. 2006 J. Chem. Phys. 125 034503 -
Vega C. , de Miguel E. 2007 J. Chem. Phys. 126 154707 -
Pi H. L. , Aragones J. L. , Vega C. , Noya E. G. , Abascal J. L. F. , Gonzalez M. A. , McBride C. 2009 Mol. Phys. 107 365 -
Vega C. , Abascal J. L. F. , Conde M. M. , Aragones J. L. 2009 Faraday Discuss 141 251 -
Harris K. R. , Woolf L. A. 2004 J. Chem. Eng. Data 48 1064 -
NIST Chemistry WebBook
Wu Y. J. , Tepper H. L. , Voth G. A. 2006 J. Chem. Phys. 124 024503 -
González M. A. , Abascal J. L. F. 2010 J. Chem. Phys. 132 096101 -
Wensink E. , Hoffmann A. , van Maaren P. , van der Spoel D. 2003 J. Chem. Phys. 119 7308 -
Hess B. 2002 J. Chem. Phys. 116 209 -
Zhang Y.-G. 2001 Mol. Phys. 99 283 -
Chen T. , Smit B. , Bell A. T. 2009 J. Chem. Phys. 131 246101 -
Tazi S. , Boþan A. , Salanne M. , Marry V. , Turq P. , Rotenberg B. 2012 J. Phys.: Condens. Matter 24 284117 -
Nosé S. 1984 J. Chem. Phys. 81 511 -
Hoover W. G. 1985 Phys. Rev. A 31 1695 -
Swope W. C. , Andersen H. C. , Beren P. H. , Wilson K. R. 1982 J. Chem. Phys. 76 637 -
Evans D. 1977 J. Mol. Phys. 34 317 -
Evans D. J. , Murad S. 1977 Mol. Phys. 34 327 -
Lee S. H. 2013 Bull. Korean Chem. Soc. 34 3800 -
Lee S. H. , Park D. K. , Kang D. B. 2003 Bull. Korean Chem. Soc. 24 178 -
Lee S. H. 2004 Bull. Korean Chem. Soc. 25 737 -