The Radial Distribution Functions of the Scaled OSS2 Water

Abstract

Classical molecular dynamics (MD) simulations using a scaled OSS2 potential originally derived from ab initio calculations are used to study the radial distribution functions of water. The original OSS2 water potential is shown to represent a glassy or an ice at ambient temperature, but the diffusion coefficient increases on increasing the temperature of the system or decreasing the density. This suggests scaling the OSS2 potential. The O-O, O-H, and H-H radial distribution functions and the corresponding coordination numbers for the scaled OSS2 potential, obtained by MD simulation, are in good agreement with the experiment results and calculations for the SPC/E water potential over a range of temperatures.

INTRODUCTION

In the study for the dynamics of H+ in water, a dissociable water potential is essential to describe how water solvent molecules can participate in ionic chemistry through dissociation and reassociation of H+ in OH−, H2O, and H3O+. Several attempts at using dissociating water potentials in classical simulations have been made previously beginning with pioneering work on central force potentials by Stillinger, Lemberg, David and Weber.1-5 In their potential the water molecules consist of H+ and polarizable O2− ions whose bare Coulomb interactions are modified at short range.

Recently Ojäme et al.6 reported progress in the design of a family of potentials for describing H+(H2O)n, called OSS (Ojäme-Shavitt-Singer)n (n=1−3). The potentials were generated by fitting to results of ab initio electronic structure calculations for the H5O2 + ion, the H2O molecule, and the H3O+ ion, as well as some results for the neutral water dimer to an analytic pair-potential. The potential could well reproduce ab initio results for the H5O2 + ion, and also provide formation energies and structures of both protonated- water and water-only clusters that agree favorably with ab initio Møller Plesset (MP2) calculations.

We have chosen the OSS2 potential for the previous studies since it was a preferred choice for MD simulation studies, because of the faster and less elaborate computercode implementation as compared to the OSS3 potential, even though it is reported that the best results were obtained using the OSS3 potential and that the OSS2 potential also gave good results, but usually exhibited too large bond angles for water molecule.6 We have reported some results7,8 using the OSS2 potential. We found, however, that the original OSS2 potential behaves like a glass or an ice under ambient conditions8, as discussed with L. Ojamäe9 and also suggested by N. Agmon.10 We also found that glass formation is absent at higher temperature (~500 K) or lower density (~0.5 gcm−3).8

In the present paper, the main goal is to seek a relevant way to eliminate the glass formation of the OSS2 water potential at ambient conditions. We accordingly scale the potential to agree with the diffusion coefficient under ambient conditions and compare the atom–atom radial distribution functions of the scaled OSS2 potential with the experimental results11 and results from simulation studies of other non-dissociating potentials for water.

Molecular Potentials and MD Simulation Methods

In the OSS2 potential potential, the total energy is given by

The first term represents the total electrostatic energy,

where rij=ri−rj, and Tij is the dipole tensor. Here nO and nH are the number of oxygen and hydrogen atoms, respectively, qi is the charge on particle i (+e for hydrogen and −2e for oxygen), μi is the induced dipole on oxygen i and α is its polarizability, and are the electric field cutoff functions for charge-dipole and dipole-dipole interactions, respectively. The induced dipole moment at each oxygen site can be obtained self-consistently by imposing the conditions dVel/dμk=0, k=1,2,…, nO:

The second and third terms of Eq. (1) represent pairwise additive potential-energies between the H and O atoms and between the O and O atoms, respectively. In addition to the electrostatic and pairwise additive terms, the last term of Eq. (1) represents a three-body term. This term is short range and describes the interaction within H−O−H triplets.

We used Gaussian isokinetics12-15 to keep the temperature of the system constant. Ewald summations were used in our simulations with the parameter for κ=5.0/L and the real-space cut distance rcut and Kmax chosen as 0.5L and 7, respectively, where L is the length of the box. (L=1.864 nm for 216 water molecules of this study) The double summations in reciprocal space, which cannot be reduced to a single summation due to the cutoff functions, were ignored. This is reasonable as the distances in reciprocal space are larger than the length L of the box. The equations of motion were solved using velocity Verlet algorithm16 with a time step of 10−15 second (1 fs). The simulations were first validated by checking our results against Ojamäe’s work for pure water using the same OSS2 potential. The calculated oxygen-hydrogen (O−H) radial distribution function g(r) and the hydration number n(r) for hydrogen in the 216 molecule pure water system are nearly identical,7(a) even though Ojamäe et al. used a different method6 for the Ewald summations in the calculation of the induced dipole moment. The equilibrium properties at each temperature are averaged over five blocks of 2,000,000 time steps, for a total of 10,000,000 time steps (10 ns) after for 2,000,000 time steps to reach an equilibrium state. The configuration of each ion is stored every 10 time steps for further analyses.

Scaling and Results

The OSS2 potential was developed from ab initio calculations of proton transfer between pairs of water molecules and it is not surprising that collective many body effects of water molecule in the condensed phase would require modification of the pair potential in condensed phase. Since the potential is multiplied by the inverse temperature in the exponential term of the partition function, scaling the potential at constant temperature (T) is equivalent to multiplying the inverse temperature by the same factor at constant potential. We found empirically that the diffusion coefficient of water at a given high temperature (T') using the original potential is the same as the experimental diffusion coefficient of water at room temperature. We have accordingly scaled the potential to agree with the diffusion coefficient under ambient conditions hoping to capture the many-body effects by this effective pair potential.

The scaled OSS2 potential is defined by VsOSS2 = λ VOSS2, where the parameter λ determines VsOSS2 from VOSS2 which is given as Vtot in Eq. (1). This requires that the charge qi of particle i in Eq. (2) is scaled by the square root of λ (qi' = λ1/2qi) in calculations of the electrostatic energy, while the induced dipole moment at each oxygen site is obtained self-consistently in the same way as it was before scaling.6,8 The other three non-electrostatic potential energies (VOH, VOO, and VHOH) in the OSS2 potential are also scaled by λ.

The diffusion coefficient for the OSS2 potential at room temperature and fluid density r = 0.9970 g/cm3 is near zero8, but at 540 K and the same density the diffusion coefficient for the system of N = 216 water molecules is nearly identical to the experimental value for water (2.26- 2.29 × 10−5 cm2/sec 17-19) at 298.15 K. From our scaling hypothesis λ = T/T' = 298.15/540 = 0.552, and we expect the scaled OSS2 potential at 298.15 K to have nearly the same diffusion coefficient (2.30 ± 0.09 × 10−5 cm2/sec) as the OSS2 potential at 540 K. The diffusion coefficient with λ = 0.552 at 298.15 K from the mean square displacement (MSD) (2.00 ± 0.09 × 10−5 cm2/sec) is within 12% of the experimental result. Further fine-tuning by choosing λ = 0.530 leads to agreement (D = 2.27 ± 0.07 × 10−5 cm2/sec) to within 1%.20 For the system of N = 32 water molecules, we found that λ = 0.801.21

The atom−atom O−O, O−H, and H−H distribution functions for the scaled OSS2 potential (VsOSS2=0.530 VOSS2) are shown in Figs. 1(a), 1(b), and 1(c), which are compared with those obtained from the original OSS2 potential, the experimental results for liquid water from neutron and X-ray diffraction data11, and distribution functions for the rigid SPC/E potential for water22 at 298.15 K. The peak heights are lowered and the valley heights are raised in moving from the OSS2 to the scaled OSS2 potential at room temperature but the positions are nearly the same. The distribution functions of the sOSS2 potential are in good agreement with experiment in those three distribution functions.

Fig. 1.(a) O−O, (b) O−H, and (c) H−H radial distribution functions at 298 K. Solid line: the experimental result, dotted line: the SPC/E potential, dashed line: the original OSS2 potential, and long-dashed line: the scaled OSS2 potential with λ=0.530 for the system of N=216 water molecules.

Fig. 2.(a) O−O, (b) O−H, and (c) H−H coordination numbers at 298 K. The legends are same as in Fig. 1.

When compared with those from the SPC/E potential, the results for the scaled OSS2 potential are comparable or superior: for the gOO(r) in Fig. 1(a), the distribution functions are located in the opposite directions each other from the experimental result with less deviation of the SPC/E result. In Figs. 1(b) and 1(c), the gOH(r) and gHH(r) of the SPC/E potential at the first peak are a kind of delta function because of the rigidity of OH bond. The gOH(r) of the scaled OSS2 potential are superior to those of SPC/E potential at the second peak, while it is opposite for the gOH(r) at the third peak and for the gHH(r) at the second and third peaks.

Figs. 2(a), 2(b), and 2(c) show the atom−atom O−O, O−H, and H−H coordination numbers of the scaled OSS2 potential, which are also compared with those obtained from the original OSS2 potential, the experimental results for liquid water from neutron and X-ray diffraction data,11 and those for the rigid SPC/E potential for water22 at 298.15 K. The nOO(r), nOH(r), and nHH(r) of the scaled OSS2 potential are much improved from those of the original OSS2 potential. The results for the scaled OSS2 potential are also comparable or superior to the SPC/E potential with step functions for the nOH(r) and nHH(r) of the SPC/E potential due to the rigidity of OH bond.

Finally we compare the atom−atom O−O, O−H, and H−H distribution functions for the scaled OSS2 potential at 268 K and at 423 K with the experimental results for liquid water from neutron and X-ray diffraction data11 in Figs. 3(a) and 3(b). Generally speaking, the agreement is satisfied except some points such as the first peaks of gHH(r) at 268 K and gOO(r) at both temperatures, and the second valleys of gHH(r) at both temperatures. The agreement of gOH(r) at both temperatures is rather considerably good over the whole range.

Fig. 3.(a) Radial distribution functions at 268 K and (b) at 423 K. Solid lines: the experimental results and the other lines: the scaled OSS2 potential with λ=0.530 for the system of N=216 water molecules.

 

CONCLUSION

Classical molecular dynamics (MD) simulations using the OSS2 potential derived from ab initio calculations can be used to understand the dynamic behavior of the OSS2 water at several different states. We found that reducing the total potential energy as well as decreasing the density or increasing the temperature of the system invokes the vivid mobility of the OSS2 water molecules. By choosing the scaling factor λ of the total potential energy as 0.530 for 216 water or equivalently scaling the charge of particle by λ1/2 = 0.728, the obtained MD results for the radial distribution functions and the corresponding coordination numbers of the scaled OSS2 water at several different temperatures are in general agreement with the experimental results, which indicates the original OSS2 potential energy is too strong resulting in the behavior like a glass or an ice.