Molecular Dynamics Simulation Study for Shear Viscosity of Water at High Temperatures using SPC/E Water Model

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 rcut and Kmax chosen as 10.0 Å and 7, respectively. Nosé-Hoover thermostat30,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 algorithm32 with a time step of l0-15 second (1 fs) and a quaternion formu-lation33,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).

Table 1.aRef.22

The original GK formula for viscosity ƞ is given by

where

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.

Figure 1.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.

Figure 2.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 = Do exp(−ED/RT)] are 3.1 and 3.7 kcal/mol for DMD and for DExp, respectively.35

According to the Stokes-Einstein (SE) relation about a Brownian (B) particle immersed in a solvent (s), DB = kBT/ CπƞsrB, where C is the hydrodynamic boundary condition, and ƞs and rB are the viscosity of solvent and the radius of the B particle. Applying the SE relation to bulk water, ƞD ∝ kBT. 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/200516 at high temper-atures are currently under investigation.

Figure 3.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 K35 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.