Temperature Dependence on Structure and Self-Diffusion of Water: A Molecular Dynamics Simulation Study using SPC/E Model

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.

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) (original1 and modified2), SPC (simple point charge) (original3 and refined4), 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 conductivity31-34 using various water models have been reported. Our interest in this study concentrated on the radial distribution functions, gOO(r), gOH(r), and gHH(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 m2/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 m2/s),16 for the modified TIP3P water model between 2.3 and 5.2 (× 10−9 m2/s),17-21 for the original SPC water model between 3.6 and 5.2 (× 10−9 m2/s),16 for the refined SPC water model between 4.2 and 4.4 (× 10−9 m2/ s),22 and for SPC/E water model between 2.2 and 4.4 (× 10−9 m2/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 water3940 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) model5 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

where σ = 3.169 Å and ε = 0.6502 kJ/mol are Lennard-Jones (LJ) parameters between oxygen atoms on different water molecules, qi is the charge at site i in water (qO = −0.8476 e and qH = 0.4238 e), and rij 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 rcut and Kmax chosen as 10.0 Å and 7, respectively. Nose-Hoover thermostat4142 was used to keep the temperature constant (the Nose-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-directions, and the minimum image convention for pair potential were applied. The equations of translational motion were solved using a velocity Verlet algorithm43 with a time step of 10−15 second (1 fs) and a quaternion formulation4445 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).

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

and the Green-Kubo relation from the velocity autocorrelation (VAC) function.46

 

Results and Discussion

The radial distribution functions, gOO(r), gOH(r), and gHH(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 data47 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 gOO and gOH 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 gOO(r) and gOH(r) are shifted to longer distances with increasing temperature except those of the first minima of gOH(r) which are never changed with temperature, as seen in Table 2.

Table 2.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 gOO(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 gOO(r). For gOH(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 gOH(r). The agreement between MD result of gHH(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.

Figure 1.MD results of gOO(r) at 300-500 K and the experiment values at 298 and 423 K.

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

Figure 3.MD results of gHH(r) at 300-500 K and the experiment values at 298 and 423 K.

Figure 4.Mean square displacements (MSD) of SPC/E water at 300-550 K.

Figure 5.Velocity auto-correlation (VAC) functions of SPC/E water at 300-550 K.

Two sets of experimental measures374849 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 form37:

where a0 = 6.11903, a1 = −1.195593, a2 = −0.05229, and a3 = −0.01841 for D(Exp1) over 278-498 K and a0 = 13.2172, a1 = −9.08602, a2 = 2.80883, and a3 = −0.35713 for D(Exp2) over 242-498 K.

Figure 6.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 K38 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 study23 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:

where D0 is the pre-exponential factor, RT has the usual meaning, and ED 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 K363751 and for over the range of temperatures 243 to 298 K49 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, gOO(r), gOH(r), and gHH(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.