Interstitial Voids and Resultant Density of Liquid Water: A First-Principles Molecular Dynamics Study.

Many anomalous properties of water can be explained on the basis of the coexistence of more than one density states: high-density water (HDW) and low-density water (LDW). We investigated these two phases of water molecules through first-principles molecular dynamics simulations using density functional theory (DFT) in conjunction with various van der Waals-corrected exchange and correlation functionals. Different density regions were found to exist due to the difference in short-range and long-range forces present in DFT potentials. These density regions were identified and analyzed on the basis of the distribution of molecules and voids present. We defined a local structure index to distinguish and find the probability of occurrence of these different states. HDW and LDW arise due to the presence of "interstitial water" molecules in between the first and second coordination shells. The population of interstitial water molecules is found to affect the overall dynamics of the system as they change the hydrogen bond pattern.

Introduction

Molecular dynamics simulation studies have obtained great attention in predicting and understanding the properties of molecules through various analytical and computational methods.[1] In earlier days, classical molecular dynamics (CLMD) simulations, dependent on the development of suitable force fields, were used to interpret accurate results as compared to data obtained from experiments. In the last few years, first-principles molecular dynamics (FPMD) simulations have been explored as a tool for assessing structural and dynamical properties of the systems. The results are encouraging and have been shown to provide new opportunities though they are still limited to relatively short time scales due to the requirement of enormous computational time as compared to that in classical molecular dynamics. Kohn–Sham (KS) formalism[2,3] of density functional theory (DFT) has emerged as the preferred electronic structure method of choice for water to biomolecules because of its favorable performance-to-cost ratio, but the performance of DFT with generalized gradient approximation (GGA) exchange/correlation functionals on liquid water is not satisfactory. GGA functionals make liquid state of water overstructured and less dense, and water diffuses too slowly compared to that in experimental results.[4] Many approximate exchange/correlation functionals have failed to describe the fine balance between relatively strong directional H-bonding and nondirectional van der Waals (vdW) interactions of liquid water. The poor description of the dispersion forces in the GGA functionals is believed to be one of the reasons behind the origin of such differences. Most of the approaches currently followed to add the dispersion corrections to existing DFT functionals include various empirical components: van der Waals density functional,[5] DFT-D methods,[6−8] and dispersion correcting atom-centered one-electron potentials (DCACP).[9] Among these methods, the last two have been able to predict the density of water molecules satisfactorily. The structure of the hydrogen bond network in water is found to be affected by the type of van der Waals (vdW) correction and GGA exchange/correlation terms used, as the hydrogen bonds are formed due to the delicate balance between covalent and electrostatic effects.[1,2] The difference in the change of strength of the hydrogen bonding interaction leads to different local tetrahedral arrangements of the molecules giving rise to differences in the density, structure, and distribution of voids of the system, which further result in the difference in dynamic properties. In fact, water can be considered a “mixture” of two or more binding states,[10−12] namely, high-density and low-density water (LDW) molecules. High-density water (HDW) and LDW are so named because they give rise to different density regions in the liquid state and consequently have different hydrogen bond patterns.[13−16] The different types of water molecules arise due to the presence of interstitial molecules in extended network of the local tetrahedral arrangement. The formation of interstitials is due to both distortion and collapse of first and second solvation shells, respectively.[17] Two states of water molecules are very common and found to exist in the supercooled region of water[18−21] and also when ambient and supercooled water undergo transition under compression. Many anomalous properties of water, such as exhibition of maximum density at 4 °C, presence of a second critical point at low-temperature water, and stability-limit conjecture,[22−26] which make it so unique liquid, can be explained on the basis of the coexistence of these two phases.[27,28] Experimentally, the signatures of these different phases of water have been observed in the neutron scattering experiments,[12,28] terahertz spectroscopy,[27] and ultrafast spectroscopy.[11,29] The involvement of the different van der Waal’s potentials can lead to different populations of these two states.[10,30,31] Therefore, it would be interesting to observe the effect of these different model potentials having different van der Waals and GGA exchange/correlation terms on the structure and density, distribution of interstitial molecules, and distribution of void space of liquid water.

In this respect, we performed FPMD simulations having different van der Waals and GGA exchange/correlation terms to probe the liquid phases of water. We considered Becke–Lee–Yang–Parr (BLYP), BLYP-D2, BLYP-D3, BLYP-DCACP, Perdew–Burke–Ernzerhof (PBE), PBE-D2, PBE-D3, and PBE-DCACP model potentials to simulate water at 300 K. All simulations were carried out in the isobaric–isothermal ensemble to allow the systems to find their equilibrium densities. Sampling within the isobaric–isothermal ensemble has proved to be an invaluable tool for accessing the relatively accurate interaction potential. This allows one to directly reproduce the experimental densities under various thermodynamic conditions by using an appropriate combination of electronic parameters. The convergence of pressure and the low-frequency volume fluctuation are the two key reasons for not adopting this method for simulations of a larger number of particles. Our main aim is to determine how these model potentials affect the distribution of HDW and LDW populations and whether that can be used to describe the dynamics of the system. We identified the HDW and LDW on the basis of the radial distribution function (RDF), spatial distribution function (SDF), and local structure index (LSI). We introduced a local structure index parameter to calculate the inhomogeneity near the end of the first solvation shell and beginning of the second solvation shell to give us an overview of the presence of interstitial water molecules. Calculation of the voids actually provides us with the structural details of a system and the information about the extended hydrogen bond network. It is well known that the environment in which a water molecule is residing can be well described by the vibrational spectral diffusion of hydroxyl groups of the water molecules.[32] This is because the vibrational frequencies of these modes are directly related to the strength of the bond, which is again dependent on the surrounding of the water molecules. Therefore, the results of spectral diffusion can be used to correlate our analysis of HDW and LDW populations with the overall dynamics of the system. Our work is important in the sense that the effect of first-principles model potentials on the hydrogen bond pattern in liquid water was analyzed in terms of the presence of HDW and LDW, which can explain the dynamical properties.

Results and Discussion

We first explored the liquid density of water with various electronic parameters from the initial isothermal–isobaric simulations. With the addition of dispersion correction, the density of liquid water enhanced drastically. The density of BLYP-D2 is very close to the liquid density of water. We obtain densities 765, 1021, 1095, and 1082 kg m–3 for simulations with BLYP, BLYP-D2, BLYP-D3, and BLYP-DCACP, respectively. Overall, BLYP functional provided an underestimated density than that of PBE functional. The densities 792, 1026, 1108, and 1091 kg m–3 were obtained for simulations with PBE, PBE-D2, PBE-D3, and PBE-DCACP, respectively. The obtained densities from these methods show a remarkable spread with respect to the experimental value due to improper description of interactions between water molecules by functionals without dispersion corrections. The corrected interactions provide cohesive interactions between water molecules in the liquid state. As found for liquid water from earlier studies,[4,33,34] the density of liquid water is quite sensitive to the exchange/correlation functional, and the noncorrected functionals estimate densities, which are around 20% lower than the experimental one. We note that the water density from the strongly constrained and appropriately normed density functional based ab initio theoretical study is ≃5% larger than that of the experimental value, and this value is a significant improvement over the PBE and BLYP functionals.[35] The analysis of structural arrangement of water molecules with various dispersion corrections was investigated by calculating the oxygen–oxygen radial distribution function (RDF) (Figure ). The characteristics peak features depend upon the choice of different functionals and the choice of the dispersion corrections. The BLYP, BLYP-D2, BLYP-DCACP, PBE-D2, and PBE-DCACP functionals overestimate the oxygen–oxygen RDFs. These functionals showed higher intensity and stronger first maxima, followed by a much deeper minima position compared to the experimental data. The positions corresponding to first maxima of oxygen–oxygen RDFs of PBE and PBE-D3 are very close to those of the experimental data, whereas the first minimum position of oxygen–oxygen RDF is in good agreement for BLYP-D3. However, the peak height is shorter than that in the experimental data.[36,37] For better understanding the RDFs, we have also calculated the spatial distribution functions (SDFs) in which the oxygen clouds are shown around the water molecule in the first and second solvation shells. The calculated results for SDFs are shown in Figure . It is found that there is a gap between first and second solvation shells of oxygen cloud density. The cloud density of oxygen around the water molecule in the first solvation shell is much denser in the case of dispersion corrected PBE functionals with respect to the BLYP, which may be the reason for the slower dynamics liquid water having former functionals.[38,39] On further observation, it is seen that BLYP forms a clear two-coordination shell, whereas its dispersion corrected functionals do not have a well-formed second coordination shell, and they are more distorted or diffused. In case of the PBE and its dispersion corrected functionals, PBE-DCACP seems to have little distorted two-coordination shells. The presence of the diffused and distorted second coordination shell is an indication of the presence of some water molecules, which are not well placed in the first or second coordination shell and are present in between the two coordination shells.[40]

Figure 1

Radial distribution functions (solid lines) for oxygen–oxygen pairs for BLYP and PBE functionals are shown in upper and lower panels, respectively. The experimental data for oxygen–oxygen RDF are shown from neutron scattering experiment by Soper at 298 K.[36]

Figure 2

Spatial distribution functions (SDFs) of oxygen cloud around water molecules. The red and yellow colors represent the oxygen cloud in the first and beyond the first solvation shell, respectively.

The gap observed in distribution of oxygen clouds corresponding to the first and second solvation shells of SDFs can give us an overview of the structurally different water molecules, which primarily is the reason for the existence of HDW and LDW. The main structural difference between these two categories is due to the presence of the interstitial water molecules, which can be better understood if we plot the higher order distributions.[28,41,42] Interstitial water molecules are the water molecules that are mainly trapped between the first and second solvation shells. The quantification of these water molecules can be done by many ways. Any change in temperature and pressure or in the model potential will affect the hydrogen bond interaction[30] between the water molecules and will shift the equilibrium between these two forms in a continuous way. In HDW, the first tetrahedral shell remains unaffected, but the second shell approaches the interstitial sites of the first shell, resulting in HDW structures. There is a possibility of significant overlap of the first and second solvation shells; more precisely, it can be said that the molecules labeled as instantaneous first shell neighbors may be second shell neighbors. Under ambient condition, the ∠OI–O–OI distribution, where OI is the oxygen molecule present at the first solvation shell of the central oxygen O, peaks at about 104.5°. Ideally, the ∠OI–O–OII angle distribution, OII being the oxygen molecule present in the secondary shell of O, gives two peaks at 45 and 75° under ambient pressure, which merges to 65° with the increase in HDW population. A representative model regarding the position of water molecules with the angle distribution is presented in Figure . This shift is due to the change in the structural topology of the second solvation shell, where the other second shell molecules are closer to the central molecules than they would have been under ambient pressure. High-density second shell molecules are not hydrogen bonded; they are purely interstitial nonbonded water molecules, and generally they make angles closer to 65°. To structure out the distribution of the interstitial molecules, we plotted the angle distribution of ∠O–OI–OII with respect to the O–OI and OI–OII distances involved in Figure . Generally, these water molecules are found in the range between 2.7 and 3.35 Å, that is, the gap between the first and solvation shells.[13] It varies from model to model we are using, depending upon the first minima in the RDF. The more well defined the RDF is, the less probability of finding interstitial water molecules. The presence of interstitial water molecules will decrease the angle between the neighbors.[40] Therefore, a careful observation of the ∠O–OI–OII in the range of 2.7–3.35 Å (i.e., the region where the interstitial molecules are generally present) can give us an idea of the distribution of the interstitial water molecules present in the system. In Figure , we have plotted the distribution of ∠O–OI–OII with the corresponding O–OI and OI–OII distances. The x axis and the y axis show the O–OI and OI–OII distances, respectively, and the color bar legend corresponds to the calculated angle between O–OI–OII. Moreover, the different colors of the contour plot denote the probability of finding a particular value of ∠O–OI–OII. It can be seen that at higher oxygen–oxygen distances, in the range of 2.7–3.35 Å in the x and y axes, where the probability of finding interstitial molecules is more, the top right-side corner of the contour plot is blue for all of the functionals. This blue color corresponds to probability of finding lower angles. This means that as the oxygen–oxygen distance is increasing (toward the end of the first solvation shell), the probability of finding some interstitial water molecules, which make the angles smaller than the tetrahedral angle, is more. When compared with different water models, BLYP is found to have a lesser blue region than that of the dispersion corrected BLYP functionals. This indicates lesser population of HDW in BLYP; however, we found PBE-DCACP to have a considerable blue region at this distance. Therefore, both the choice of the exchange/correlation functional and the dispersion corrections are very sensitive toward the population of both types of water densities.

Figure 3

Representative model of the first solvation structure of water (O and OI atoms) with few of the molecules from the second solvation shell (OII atoms). The green and brown ones depict the hydrogen-bonded and interstitial water molecules, respectively.

Figure 4

Contour plots of O–OI and OI–OII distances involved with ∠O–OI–OII, as shown in Figure . The x axis and y axis represent O–OI and OI–OII distances, respectively. The color ranges given in the color bar legend are according to the values of angles.

To have a better understanding about the populations of the HDW and LDW water, we calculated a local structure index, S, which is defined aswhere S(i,t), the local structure index of ith oxygen at time t, shows the distribution of oxygen atoms placed in the second solvation shell of the ith oxygen with respect to the first solvation shell. Therefore, d(i,j) is the distance between the ith oxygen atom with its second shell jth oxygen atom neighbor and rmin is the average value of the first solvation shell. n(i,t) is the number of jth neighboring atom of ith oxygen at time t. The first minimum of the g(r) corresponding to the first solvation shell can give an idea of rmin if the first and second solvation shells are well defined. In Figure , we have plotted the probability distribution, P(S), with the distance for different water models. A low value of S(i,t) will indicate a disordered local environment, that is, the presence of interstitial water near the first solvation shell whose distance is less than that of the second solvation shell. A high value of S(i,t), on the contrary, would indicate an ordered system, which means the second neighboring water molecules are more likely to be present at where the second hydration shell is supposed to be. Therefore, a careful analysis of the probability distribution function, P(S), can give us an idea of the presence of HDW population and LDW water in different water models. Figures and 5 show that apart from BLYP, all three models, BLYP-D2, BLYP-D3, and BLYP-DCACP, have considerable probability of finding water molecules at the described distance distribution, which indicates that they have interstitial water in between the second and first solvation shells. Among these three models, BLYP-D3 has the highest probability and BLYP has the least (Figure ). A careful investigation of the RDF for BLYP-D3 shows that there is a considerable flat region after the first coordination shell. In case of PBE and its dispersion corrected models, PBE-DCACP, PBE-D3, and PBE-D2 show considerable probability of finding interstitial water molecules with respect to PBE. These results are well correlated with the SDF presented in Figure . Because BLYP has well-defined first and second solvation shells, there is less probability of finding interstitial water molecules; thus, less blue region is observed in the right-hand top corner of the contour plot (Figure ). For PBE-DCACP, the SDF is visibly distorted, showing considerable presence of interstitial water.

Figure 5

Probability distribution of local structure index (S) of the secondary solvation shell oxygen molecules with respect to the first minimum of oxygen–oxygen RDF (first solvation shell).

The presence of interstitial water molecules also affects the interatomic void space by modifying the network structure. Here, we analyzed the voids in water molecules through the Voronoi polyhedra (VP) method, which is one of the efficient ways by which we can analyze the structural details of a system. Specially, it will be important for the current study, where we are comparing different models of water differing in the dispersion energy terms. The arrangement of the water molecules in the coordination sphere differs a lot depending on the model potential. From a set of configuration of atoms, we can generate a tessellation in the three-dimensional space, known as Voronoi tessellation, named after the scientist who invented the method.[43] The tessellation consists of repetitive units of VP region of the atoms. The units are defined as the regions consisting of all points that are closer to the reference atom than to any other atoms. The void radius is defined as the distance between the vertex of the VP and the atom minus the atom radius. The other way of defining the tessellation is Delaunay tessellation. VP and Delaunay tessellation are complementary to each other. We considered the coordinates of the oxygen atoms to generate the tessellation, and the Voronoi analysis was executed by the algorithm described in refs (43−47). The plot of the probability distribution of the voids present in the different systems is presented in Figure . It is seen that the methods have quite a large impact on the void distribution of the system. This can be actually related with the RDF of the system. Calculations of the distribution of the voids and radial distribution function, which actually calculates the distribution of neighboring atoms around a fixed atom, are complementary to each other. Figure shows that BLYP has more void radii than BLYP-D2 > BLYP-D3 > BLYP-DCACP. This implies that at a fixed distance, the water molecules in the BLYP model will be less surrounded by neighboring water molecules than will the BLYP-D2 model, which is further lesser than BLYP-D3 < BLYP-DCACP. This is also reflected in the calculation of the coordination sphere of the different models. In other words, we can say that BLYP water models are less dense than the other BLYP dispersion corrected models. It can be also noted here that on the basis of our Local structure index calculation, BLYP has more probability of having LDW molecules than that of the other models. Similar types of comparison can be made in cases of PBE, PBE-D2, PBE-D3, and PBE-DCACP. It is found that PBE-D2 has the least void radius, followed by PBE-DCACP < PBE ∼ PBE-D3. Further, RDF and coordination number also reflect the same scenario. Two water models, PBE and PBE-D3, have quite close coordination numbers, and the other two water models, PBE-DCACP and PBE-D2, have close coordination numbers that correlate well with the probability distribution of the void radius. Therefore, it can be said that PBE and PBE-D3 are less crowded water models than are PBE-DCACP and PBE-D2 in the first coordination shell. PBE-DCACP has more HDW population than PBE. Further, we matched our result with the spectral diffusion of the water molecules. It is found that vibrational frequency of a water molecule is sensitive to its local solvation environment. We correlated the data reported by Klein and co-workers[39] for further discussion. The dynamics of spectral diffusion can be calculated by frequency–time correlation and spectral hole dynamics calculations.[32] Both of these functions have a short decay time scale (∼100 fs) and a long decay time scale (∼2 ps) corresponding to the hydrogen bond lifetime. The slower time scale can be related to the dynamics of the local structural relaxation. The dynamics of the hydrogen bond making is found to have a faster time scale of ∼100 fs; hence, this can be correlated to the short-time dynamics of the spectral diffusion. The long-time dynamics can be related with the hydrogen bond dynamics of the system. On carefully comparing the three systems, BLYP, BLYP-D2, and BLYP-D3, it can be seen that the τ2 is decreasing in the order BLYP > BLYP-D2 > BLYP-D3. These can be explained in terms of the interstitial water that BLYP has the least population of interstitial water and BLYP-D3 has more, so the hydrogen bond making and breaking become more feasible for BLYP-D3, and therefore, the lifetime decreases. For the same reason, BLYP-DCACP shows faster dynamics than that of BLYP and PBE shows slower dynamics than that of PBE-DCACP.[38]

Figure 6

Probability distribution of radius of vacancies of water molecules for different model potentials.

Conclusions

To summarize, we have explored the effects of dispersion interactions on the structure, density, voids, and interstitial molecules of liquid water from the analysis of trajectories generated from first-principles molecular dynamics simulation using CP2K at room temperature. We presented the structural details of liquid water in terms of HDW and LDW for BLYP as well PBE and their dispersion corrected functionals. Addition of dispersion corrections to these functionals produced more dense liquid as the oxygen–oxygen first coordination number is increased, which leads to small void spaces. However, the first and second solvation shells are not well defined on addition of dispersion correction in case of BLYP functionals. These diffused coordination shells result due to the presence of interstitial water molecules between the first and secondary shells. We have plotted the spatial distribution function (SDF), where it is seen that BLYP has more compact coordination shells compared to that of BLYP-D3 and PBE has more compact secondary shells compared to that of PBE-DCACP. The presence of these interstitial water molecules gives rise to high-density water in the model functionals, and they are found to affect the dynamics of the system as they help in breaking and making of the hydrogen bonds quickly. We have calculated the population of the high- and low-density water molecules for the different functionals. BLYP-D3 was found to have the maximum population of high-density water molecules compared to that of other BLYP functionals, and PBE-DCACP was found to have more HDW water population than that of PBE. It was also concluded that PBE shows slower dynamics than that of BLYP functionals. It was seen that we can describe the time fluctuation correlation of a hydroxyl bond in the presence of HDW and LDW populations. In this study, our effort is to provide a quantitative overview of the presence of interstitial water molecules in different model potentials in terms of two different density states. We emphasize that our results are based on the analysis within the length scale of the second solvation shell. These observations will be different with the variation of thermodynamic conditions; we plan to address the subject in near future. Moreover, we want to note that water provides an intense challenge to both the molecular simulation and experimental techniques. A combination of accurate computational approaches, long-enough simulations, large-enough systems, and efficient analysis of trajectory data is required to represent water within its delicately balanced interactions. It is clear that these requirements within the first-principles molecular dynamics approach will be a challenge for addressing two density states of water molecules in a relatively long length scale.

Computational Details

FPMD simulations were carried out with the Quickstep[48] routine of computer simulation package CP2K (http://cp2k.berlios.de). The method developed in the code employs the Gaussian plane wave method[49] with a dual basis set formalism of Gaussian-type orbitals and plane waves for the electronic density to solve the self-consistent KS equations of DFT.[3] Energies and forces were calculated through the electronic structure part of CP2K. The electronic structure was optimized at every step (0.5 fs) of the trajectory following Born–Oppenheimer dynamics, so that the system is always kept in the electronic ground state. The DFT calculations used the Becke–Lee–Yang–Parr (BLYP)[50,51] and Perdew–Burke–Ernzerhof[52] exchange/correlation functionals with the norm conserving Goedecker–Teter–Hutter pseudopotentials for the core electrons. A triple-zeta valence plus doubly polarized basis set was used to obtain a good compromise between the accuracy and computational cost. The simulations with the BLYP and PBE functionals without any dispersion corrections have been treated as reference systems. To investigate the influence of dispersion forces, we performed additional MD simulations with different types of dispersion correction methods: Grimme DFT-D2, DFT-D3, and DCACP methods. A charge density cutoff at 600 Ry for the auxiliary plane wave basis set was used. The system with 100 water molecules was pre-equilibrated (3 ns simulation time) using SPC/E[53] forcefield in the NVT ensemble with a density of 1000 kg m–3 with a time step of 1.0 fs to obtain the starting geometry for FPMD simulations. This pre-equilibration period in CLMD was followed by a geometry optimization in CP2K before starting the FPMD simulations. The FPMD simulations were carried out in the isobaric–isothermal (NpT) ensemble at a temperature of 300 K and an external pressure of 1 atm. A massive Nosé–Hoover chain thermostat[54] was used to control the temperature, and the pressure was controlled via the barostat suggested by Mundy and co-workers.[4] An equilibration period of 10 ps was used for all systems, which was found to be reasonably long for the cell volume and system energy to stabilize. Then, each system was simulated for 25 ps to know the liquid state density. During this phase of the simulations, the volume of the simulation box was allowed to fluctuate, and the average volume was determined from the last 20 ps of the simulation. Subsequently, the simulations were carried out in a canonical (NVT) ensemble, keeping the box size fixed at the average value obtained previously for a given system. For all of the simulations, the integration time step was set to 0.5 fs. Each system was run for another 52 ps, out of which the saved trajectory at every time step of the last 50 ps was used for the calculations of various quantities.