Water under Supercritical Conditions: Hydrogen Bonds, Polarity, and Vibrational Frequency Fluctuations from Ab Initio Simulations with a Dispersion Corrected Density Functional.

We have studied the effects of dispersion interactions on the dynamics of vibrational frequency fluctuations, hydrogen bonds, and free OD modes in supercritical heavy water at three different densities by means of ab initio molecular dynamics simulations. The vibrational spectral diffusion, as described by the frequency fluctuations, is studied through calculations of frequency time correlation of stretch modes of deuterated water, and its relations to the dynamics of hydrogen bonds and free OD modes are described. In addition, some of the other dynamical, structural, and electronic properties such as diffusion, rotational relaxation, radial distribution functions, hydrogen bond and coordination numbers, and dipole moments are also investigated from the perspectives of their variation with inclusion of dispersion interactions at varying density of the solvent. Although some changes in the structural properties are found on inclusion of dispersion corrections, no significant difference in the fluctuation dynamics of OD stretching frequencies and also in other dynamical quantities of supercritical water are found because of dispersion effects. The dynamics of water molecules under supercritical conditions is very fast compared to the corresponding dynamics under ambient conditions. The large thermal effects at such a high temperature seem to take over any relatively minor changes that might be introduced by weak dispersion interaction.

Introduction

Information about the hydrogen bond dynamics of water molecules can be obtained experimentally by studying the vibrational spectral diffusion of stretch modes of water through time-dependent infrared spectroscopy.[1−16] In the liquid state, water molecules are in continuous motion that leads to the breaking and formation of hydrogen bonds, and this process also affects the vibrational frequencies of water. Various research groups have performed studies to investigate how these two phenomena are related to each other.[1−22] It is known that water molecules form three-dimensional hydrogen-bonded network in the liquid phase. In aqueous solutions, apart from water–water hydrogen bonds, water molecules also form hydrogen bonds with many other charged or polar molecules, and this hydrogen-bonded network can be perturbed by increasing the temperature of the solutions.

Water under supercritical conditions is of great importance because it behaves like a fluid of different characteristics at such high temperatures. The nature of supercritical water changes from liquidlike to a gaslike fluid on decrease of density without any phase transition. The electronic property of water changes as also the dipole moment and polarizability of water molecules when one moves from ambient to the supercritical states. The dynamical properties like hydrogen bond dynamics, rotational motion, diffusional motion of water molecules, and also the vibrational frequency fluctuations in supercritical water can be very different from those of the ambient water. To understand these dynamical properties, especially the dynamics of vibrational spectral diffusion and its relations with other relevant dynamical modes, one needs to look at the molecular dynamics of water at the supercritical state. The first experimental study of time-dependent infrared spectroscopy of supercritical water was done about a decade ago.[23,24] In that study, it was inferred that vibrational spectral diffusion in supercritical water, if any, would take place rather fast with a timescale of less than 300 fs. Theoretical studies of the hydrogen bond structure and dynamics in supercritical water were first carried out using molecular dynamics simulations with empirical potential models of water.[22,25−28] Subsequently, the connection of hydrogen bond dynamics with vibrational spectral diffusion was also investigated through ab initio simulations without involving any empirical potential models.[29] There have also been other first principles simulation studies of water which looked into various structural and dynamical behavior of supercritical water.[30,31] However, all these studies were performed within density functional theory for the electronic structure calculations using the well-known BLYP functional[32,33] without any long-range dispersion corrections. It may be noted that BLYP and other functionals belonging to the same generalized gradient approximation (GGA) category do not properly include the dispersion interactions, especially at long distances. In recent years, there have been a number of theoretical studies which have shown that inclusion of dispersion corrections[34,35] into BLYP or other GGA functionals provide better phase diagram, structure, and dynamics of water.[36−40] In view of these results, very recently, some work has also been carried out to look at the effects of dispersion interactions on the hydrogen bond dynamics and vibrational spectral diffusion of liquid water.[41−43] However, to the best of our knowledge, no theoretical study has yet been carried out on how the dispersion interactions affect the vibrational spectral diffusion and its correlations with various dynamical properties of water under supercritical conditions and also their dependence on the density of supercritical water. Here, we have performed ab initio molecular dynamics simulations of supercritical water to address these issues by using a dispersion corrected density functional.

In the present work, we have carried out ab initio molecular dynamics simulations of supercritical water by using the BLYP[32,33] functional with dispersion corrections,[34,35] which is referred to as the BLYP-D functional. Three different densities of supercritical water have been considered in the current study. The vibrational spectral diffusion has been investigated by calculating the frequency time correlation functions from the simulation trajectories. The calculation of the stretching frequencies is done by using the method of wavelet analysis while the hydrogen bond dynamics and rotational relaxation of water molecules are looked at by calculating appropriate time correlation functions. The BLYP-D functional is found to produce a higher peak in the O–O radial distribution function than the corresponding BLYP result at the supercritical state, and our results are found to be in good agreement with the results of an earlier study.[37] However, we have not observed any significant changes in the dynamics of spectral diffusion of water because the molecular dynamics of water molecules is very fast at this higher temperature compared to that of ambient water. It is found that the decay of vibrational spectral diffusion has primarily two timescales. At the higher density (1.1 g cm–3), there is a short timescale corresponding to the dynamics of free OD bonds and first inertial motion of OD bonds which is followed by a slower timescale corresponding to the hydrogen bond lifetime. At lower density, an opposite scenario is found where the longer-time component is found to correspond to the timescales of free OD modes. With inclusion of dispersion corrections, it is found that the longer timescale of vibrational spectral diffusion decreases from ∼0.2 to ∼0.1 ps at 1.1 g cm–3 and from ∼0.54 to ∼0.50 ps at 0.39 g cm–3. In addition to the dynamics of spectral diffusion, we have also looked at the average frequency–structure correlations for the higher and lower densities of supercritical water by using the BLYP-D functional.

We have organized the rest of the paper as follows. The effects of dispersion interactions on structural properties and also the relationship between the OD stretching frequency and the D···O distance are presented in section . The results of vibrational spectral diffusion of OD modes are discussed in section . The relations between the hydrogen bonds, free OD bond dynamics, and the vibrational spectral diffusion are discussed in section . The results of rotational correlation functions and translational diffusion are presented in section . The results of dipole moment calculations are discussed in section , and our conclusions are briefly summarized in section . We have the computational details of ab initio molecular dynamics simulations and frequency calculations in section .

Effects of Dispersion Interactions on the Structure, Hydrogen Bonds, and Frequency–Structure Correlations

In Figures and 2, we have shown the structural features of supercritical water for the BLYP-D functional, and the results have also been compared with those of the BLYP functional. The oxygen–oxygen radial distribution function is found to be more structured for the BLYP-D functional compared to that produced by the BLYP functional. Also, the BLYP-D results are found to be in better agreement with experimental results.[44] Our results are also in good agreement with the earlier results of Vuilleumier and co-workers.[37] It is found that the effects of dispersion interactions on radial distribution functions (O···D and O···O) become more pronounced as the density is decreased. We have calculated the hydrogen bond (HB) numbers for both functionals at the three densities. The presence of a hydrogen bond is determined by using a set of configurational criteria. If the O···D distance is less than the distance of first minimum of the intermolecular O···D radial distribution function and the O···O distance is less than the first minimum of the intermolecular O···O radial distribution functions, then a hydrogen bond exists between the two water molecules. The HB numbers are found to be 2.76 (2.89), 1.73 (1.90), and 1.25 (1.46) at densities of 1.1, 0.78, and 0.39 g cm–3 for the BLYP (BLYP-D) functional. The same for the bigger system of higher density is 2.8 (2.9) for BLYP (BLYP-D) functional.

Figure 1

Radial distribution functions and running coordination numbers for oxygen–oxygen pair correlations for supercritical heavy water at the density of 1.1, 0.78, and 0.39 g cm–3. The results of panels (a,b) are for BLYP-D and BLYP functionals, respectively. The different curves are as shown in the figures. *The experimental results[44] at the density of 0.75 g cm–3 are also shown in both the figures.

Figure 2

Radial distribution functions and running coordination numbers for oxygen–hydrogen pair correlations for supercritical water at the density of 1.1 g cm–3 (solid), 0.78 g cm–3 (dashed), and 0.39 g cm–3 (dashed–dotted curves). The results of panels (a,b) are for BLYP-D and BLYP functionals, respectively. *The experimental results[44] at the density of 0.75 g cm–3 are also shown in both the figures.

It is important to see how the frequency–structure correlation varies with the inclusion of dispersion interaction and also with the thermodynamic conditions. The contour plots of the conditional probability of finding a particular frequency for a given D···O distance are shown in Figure . The contour plots are found to be more elongated along the X-axis with decreasing densities under supercritical conditions. In all cases, substantial dispersions are found in the probability distributions which means no single instantaneous frequency could be assigned to a given D···O distance. However, the average frequency is seen to increase with increase of the D···O distance. The rate of increase is also seen to gradually decrease as the density is decreased. Also, the position of the maximum probability moves to a larger distance–higher frequency value with decrease of density, which can be linked to the weaker and less number of hydrogen bonds at the elevated temperature. The extent of correlation between the structure and OD frequency can be quantified in terms of the correlation coefficient between the OD frequency and the D···O distance. The correlation coefficient can be defined in the following waywhere x is the D···O distance denoted by R. Clearly, χR = 1.0 for a perfect linear correlation, while a zero value of the correlation coefficient χR means there is no statistical correlation between the OD frequency and the D···O distance. The values of χR are given in Table for all the systems studied here. The values of χR for supercritical water are found to be smaller than that for ambient water. Clearly, for supercritical water, the correlation between the stretch frequency and the hydrogen bond distance is far from perfect. However, the correlation coefficient is still far above zero, which means that a statistical correlation between the stretch frequency and the nearest D···O distance continues to exist to some extent even at supercritical states, especially when the density is not too low.

Figure 3

(a) Joint probability distributions of OD frequency and D···O distance for supercritical water at the density of (a) 1.1, (b) 0.78, and (c) 0.39 g cm–3. The results are for the BLYP-D functional.

Table 1

Values of the Correlation Coefficients for OD Frequency and D···O Hydrogen Bond Distance Correlations at Three Different Densities of Water at Supercritical Temperature for the BLYP-D Functionala

density (g cm–3)	temperature (K)	system size (N)	correlation coefficient (χ)
1.1	673	32	0.51 (0.58)
0.78	673	32	0.46 (0.50)
0.39	673	32	0.25 (0.33)
1.1	673	64	0.52 (0.54)
1.1	300	32	0.65 (0.70)

The values shown in brackets are for the BLYP functional.

Dynamics of Vibrational Spectral Diffusion

The dynamics of vibrational spectral diffusion of OD bonds is investigated by calculating the frequency–frequency time correlation function which is defined aswhere δω(t) is the fluctuation from the average frequency at time t. The results of Cω(t) are shown in Figure . It can be seen from Figure that the frequency correlations generally show a biphasic decay with a very fast decay within 100 fs and then a relatively slower decay extending over a longer time. The slower decay is seen to extend beyond 0.5 ps for the lower density of 0.39 g cm–3. To quantify the timescales of this biphasic decay, we have used a biexponential fitting of the calculated frequency–frequency time correlation function by the following function[29]

Figure 4

Time correlation functions of OD fluctuating frequencies of supercritical water for the density of 1.1, 0.78, and 0.39 g cm–3, respectively, for the BLYP-D functional. Results are also included for the density of 1.1 g cm–3 under ambient conditions. The different curves are as shown in the figure. The light gray curves show the biexponential fits.

The biexponential fits (eq ) to the simulation results are also shown in this figure. The corresponding timescales are included in Table . We first discuss the values of the shorter relaxation time (τ1). For the same system size of N = 32, τ1 is found to decrease from 0.10 to about 0.06 ps, whereas τ2 is found to increase from 0.13 to 0.50 ps when the density is decreased from 1.1 to 0.39 g cm–3. For the larger system of N = 64, the values of τ1 and τ2 are 0.07 and 0.16 ps, respectively. All these quoted values are for the BLYP-D functional although the BLYP values are not too different. The error bars of these timescales, which were calculated by dividing the total trajectory of a given system into blocks of 7.5 ps, were found to be less than 8% for the smaller systems of N = 32 and about 3% for the larger system of N = 64. As expected, the timescales are found to be better converged for the larger system. Combining both the density and system size effects, the slowing down of the second relaxation time (τ2) is found to be more significant than the acceleration of the short-time dynamics (τ1) with decrease of density. It is, however, noted that the weight of the short-time dynamics is significantly higher than that of the longer-time relaxation characterized by τ2. For τ1, the effect of density decrease is found to be almost as much as the system size effect, and thus, no definite conclusion could be drawn regarding the variation of this timescale. This variation is found to be as much as the scattering of its calculated values. Regarding the effect of the dispersion corrections, it is seen that the dispersion corrections make the dynamics significantly faster for ambient water, which was also found in earlier studies.[37,39−43] However, for supercritical water, the corresponding change in the dynamics due to dispersion corrections is found to be much less significant. We also note that each OD mode has been treated as a local oscillator for calculations of vibrational frequencies. Because the effects of vibrational population relaxation and Forster energy transfer through resonant excitonic couplings are not incorporated in our current calculations, the present work can be better considered as an approximation of the spectral diffusion of isotopically diluted water. In this context, we note the recent ab initio molecular dynamics simulation study[45] of vibrational spectral diffusion in liquid water including van der Waals corrections where the intermolecular energy transfer and the delocalization of the O–H stretch mode were shown to be particularly important for the vibrational spectral diffusion in neat liquid water.

Table 2

Spectral Diffusion Results for OD Modes of All the Systems Studied in This Worka

density (g cm–3)	temperature (K)	system size (N)	τ1	τ2	a1
1.1	673	32	0.10 (0.10)	0.13 (0.15)	0.90 (0.88)
0.78	673	32	0.078 (0.08)	0.18 (0.19)	0.80 (0.78)
0.39	673	32	0.058 (0.06)	0.50 (0.56)	0.69 (0.71)
1.1	673	64	0.071 (0.075)	0.16 (0.13)	0.91 (0.75)
1.1	300	32	0.08 (0.10)	1.30 (1.85)	0.66 (0.68)

The time constants (ps), frequency (cm–1), and weights of the time-dependent frequency correlations of OD bonds for the BLYP-D functional. The values shown in brackets are for the BLYP functional.

To have a better understanding about the dynamical origin of the vibrational spectral diffusion timescales, we have calculated the hydrogen bond dynamics and also the survival dynamics of free OD bonds. The rotational dynamics of OD modes are also calculated, and the results of these calculations are discussed in the next two sections.

Dynamics of Hydrogen Bonds and Free OD Modes

The calculation of hydrogen bond dynamics is done by using the so-called population correlation function approach.[46−51] In particular, we have calculated the continuous hydrogen bond correlation function, SHB(t), which gives the probability that a pair of water molecules, which was hydrogen-bonded at time t = 0, remains hydrogen-bonded continuously up to time t. The integral of SHB(t) provides the average lifetime of a hydrogen bond. We have defined hydrogen bond by using a set of configurational criteria. If the O···D distance is less than the distance of first minimum of the intermolecular O···D radial distribution function and the O···O distance is less than the first minimum of the O···O radial distribution function, then a hydrogen bond exists between the two water molecules. An OD mode is considered to be dangling or free when it is not hydrogen-bonded. Like hydrogen bond dynamics, we have also calculated the correlation function (SDB(t)) of free OD bond dynamics by using a similar population correlation function approach.[52] The results are shown in Figure a,b. The corresponding lifetimes, τHB and τDB, are given in Table . It can be seen from Table that for the higher density systems (1.1 g cm–3), τDB is less than τHB for both N = 32 and 64. Thus, for higher density, an OD bond remains free for a shorter time than it remains hydrogen-bonded even at supercritical temperature. Similar observation can also be made for ambient water (Table ) although the absolute values of the timescales are much longer for ambient water. For the lower density systems, however, the two timescales of τHB and τDB are found to follow a reverse order with τHB now less than τDB (Table ). Thus, for the lower density systems of 0.78 and 0.39 g cm–3, an OD bond remains free for a longer period compared to the period over which it remains hydrogen-bonded. Because the dynamics of hydrogen bond breaking and reformation are believed to be the primary factors which influence the OD frequency changes, the shorter timescale of the frequency fluctuations (τ1) may be linked to the shorter of the two timescales of hydrogen bond and free OD bond dynamics, that is τHB and τDB, while the longer timescale may be connected to the longer of τHB and τDB. However, a close look at the numerical values of Tables and 3 reveal that this linking of τ1 to the shorter and τ2 to the longer of τHB and τDB is not perfect. This means the dynamics of frequency fluctuations is not merely determined by the dynamics of associated hydrogen bond breaking and reformation, but other dynamical modes, possibly involving more than just the nearest neighbor, also contributing to the frequency fluctuations of an OD mode. Also, in addition to the free OD modes, the relaxation of weakly hydrogen-bonded modes is also expected to contribute to the dynamics at the higher frequency side. This is more so because of the simple geometric definition that we have used to characterize the free OD modes. A more robust geometric definition of the free O–H modes of water has been provided recently by Nagata and co-workers[53] in the context of air–water interfaces. Regarding the effects of dispersion corrections, we note that the hydrogen bond lifetime for pure water under ambient conditions for the BLYP-D functional is shorter than that produced by the BLYP functional.[54] Thus, under ambient conditions, the BLYP-D functional produces a faster hydrogen bond dynamics. For supercritical water also, the BLYP-D functional is found to produce a slightly faster dynamics compared to that for the BLYP functional, although the timescales of supercritical water are rather short, and the changes in relaxation times may be within the scattering caused by system size effects.

Figure 5

Time dependence of the (a) continuous hydrogen bond correlation and (b) dangling OD probability function for supercritical water for the BLYP-D functional. The different curves are as shown in the figures. The top solid curve in (a) is for the density of 1.1 g cm–3 under ambient conditions.

Table 3

Average Lifetimes (in ps) of Water–Water Hydrogen Bonds (HBs) and Dangling Bonds (DBs) of the Systems 1–4 for the BLYP-D Functionala

density (g cm–3)	temperature (K)	system size (N)	τHB	τDB
1.1	673	32	0.11 (0.14)	0.07 (0.08)
0.78	673	32	0.10 (0.12)	0.16 (0.15)
0.39	673	32	0.08 (0.07)	0.30 (0.40)
1.1	673	64	0.11 (0.13)	0.07 (0.074)
1.1	300	32	1.9 (2.30)	0.24 (0.12)

The values shown in brackets are for the BLYP functional.

Rotational Dynamics of OD Bonds and Translational Diffusion

We have also calculated the density dependence of the dynamics of rotational motion of OD bonds of water molecules under supercritical conditions. It is known from previous studies that rotational motion plays an important role in the breaking of hydrogen bond dynamics and hence in the spectral diffusion.[29,55−59] Therefore, it would be interesting to see to what extent the dispersion interaction affects this motion. In Figure , we have shown the decay of orientational function ⟨cos θ(t)⟩, where θ(t) is the angle between the orientation of an OD vector at times 0 and t, and the average is carried out over all the OD bonds present in a system. A biphasic decay is clearly seen in the figure. There is a fast decay due to inertial rotation of the molecules which is followed by a relatively slower decay. The timescales of the biexponential fit (τrotfast and τrotslow) and associated weights are given in Table . At low density (0.39 g cm–3), the increasing weight factor of fast inertial motion of OD mode is found to be a key factor for the hydrogen bond dynamics and spectral diffusion because the faster timescale in spectral diffusion is very much closer to the hydrogen bond lifetime and inertial rotational time at low density. At higher density (1.1 g cm–3), the hydrogen bond lifetime is found to be longer than the inertial rotational time. As the coefficient of inertial component is relatively small at higher density, the vibrational spectral diffusion of water molecules generally takes place through slower dynamics of hydrogen bond breaking. Thus, the hydrogen bond breaking is believed to be responsible for the slower component of the spectral diffusion at higher density. The inclusion of dispersion correction produces slight changes in the rotational dynamics of supercritical water. However, such changes appear to be not so significant because under supercritical conditions the molecular dynamics of water molecules is very fast due to thermal effects. In other words, the strong thermal effects take over any small changes that might be introduced by weak dispersion corrections.

Figure 6

Rotational relaxation of OD bond vectors of supercritical water for the BLYP-D functional. The different curves are as shown in the figure. The top solid curve is for the density of 1.1 g cm–3 under ambient conditions.

Table 4

Orientational Relaxation Time of OD Bonds of Water Molecules at Different Density for the BLYP-D Functionala

density (g cm–3)	temperature (K)	system size (N)	τrotfast	τrotslow	arotfast
1.1	673	32	0.09 (0.09)	0.25 (0.28)	0.33 (0.15)
0.78	673	32	0.073 (0.08)	0.20 (0.25)	0.43 (0.45)
0.39	673	32	0.069 (0.06)	0.12 (0.15)	0.72 (0.60)
1.1	673	64	0.08 (0.07)	0.30 (0.32)	0.32 (0.27)
1.1	300	32	0.09 (0.1)	8.4 (10.0)	0.08 (0.1)

The values shown in brackets are for the BLYP functional.

To see the effects of dispersion corrections on the translational motion of water molecules, we have calculated the diffusion coefficient of water molecules for both functionals by calculating the mean square displacements. The diffusion coefficient is found to be ∼30 (47.4) (×10–5 cm2 s–1) and ∼113.0 (89.27) (×10–5 cm2 s–1) for the density of 0.78 and 0.39 g cm–3, respectively, for the BLYP-D (BLYP) functional at supercritical conditions for N = 32. The current result of the diffusion coefficient at the lower density (0.39 g cm–3) can be compared with the available experimental value of 137 ± 7 (×10–5 cm2 s–1) at 0.23 g cm–3 at 673 K.[60] For the higher density system (1.1 g cm–3) with N = 64, the diffusion coefficient is found to be 20 (18.9) (×10–5 cm2 s–1). Clearly, significant dependence of the diffusion coefficient on system size is found for supercritical water when N is changed from 32 to 64. The effect of the system size on diffusion was also found to be significant for ambient water in earlier studies.[61]

Dipole Moment Distributions

At supercritical conditions, the dynamics of water molecules is very fast. The solvation properties of supercritical water also differ very much from those of ambient water. The reduction of polarity and the reduced hydrogen bond formation ability are believed to be the key factors for this altered solvation feature.[29] Thus, we have calculated the average dipole moment of water molecules for the current systems. The dipole moments of individual water molecules in the supercritical conditions at three different densities have been calculated by using the maximally localized Wannier functions and Wannier function centers method.[62,63] We have observed that the average dipole moment of a molecule decreases on reduction of density for both the functionals. Overall, only small changes in the dipole moment are found when the results are compared for the two different functionals at a particular density. The results of the dipole distributions are shown in Figure , and the values for the average dipole moments are given in Table for both the functionals. The results of dipole moment distribution of water molecules under ambient conditions have also been included in Figure for comparison.

Figure 7

Distributions of molecular dipole moments in supercritical water for the BLYP-D functional. The different curves are as shown in the figure. The magenta-solid curve at the right side is for the density of 1.1 g cm–3 under ambient conditions.

Table 5

Values of Average Dipole Moment of Water Molecules at Supercritical Temperature for the BLYP-D Functionala

density (g cm–3)	temperature (K)	system size (N)	dipole moment (debye)
1.1	673	32	2.55 (2.65)
0.78	673	32	2.30 (2.43)
0.39	673	32	2.05 (2.1)
1.1	673	64	2.52 (2.55)
1.1	300	32	2.96 (3.0)

The values shown in brackets are for the BLYP functional.

Summary and Conclusions

We have carried out first principles molecular dynamics simulations of supercritical water at three different densities using the dispersion-corrected BLYP-D functional. Calculations with the BLYP functional are also done for a system with a larger number of water molecules. It is found that at supercritical conditions, the peak height of the O···O radial distribution function of supercritical water decreases with decrease of density of pure water. The peak heights of O···O and O···D radial distribution functions for the BLYP-D functional are found to be greater than that of the BLYP functional. The difference between the peak heights of the radial distribution functions for the BLYP and BLYP-D functionals is found to increase with decrease of density. Slight changes in the dynamical properties are also found on inclusion of dispersion corrections. At supercritical temperature, the dispersion-corrected BLYP-D functional is found to produce very similar hydrogen bond dynamics as the BLYP functional for both higher and lower densities. In case of the free OD bond dynamics, the difference between the two functionals is found to be very small at higher densities, but some noticeable differences are found at the two lower densities. Also, when the system size is increased from 32 to 64 water molecules for the higher density, some noticeable changes in the structural and dynamical properties of supercritical water are found. These changes due to system size effects are generally not large but still they fall in a similar range as the changes caused by the dispersion corrections, which means further simulations of even bigger systems would be required to arrive at numbers that are converged with respect to the system size. This issue of system size effects is further discussed below. From the present BLYP and BLYP-D results for a given system size, it may be concluded that the dispersion corrections have only minor effects on the vibrational spectral diffusion and other dynamical properties of water under supercritical conditions.

We finally note that although the present study involved many ab initio simulations of supercritical water of different density, functional, and system size, still it can be further extended in many directions. The main limitations of the current study are (i) use of a rather small system size (N = 32) for most of the simulations except two systems where N = 64 was used, and (ii) use of only one dispersion corrected functional (BLYP-D) where many other choices are also available in the literature. As discussed in earlier sections, noticeable effects of system size are found in our calculations. The box length for the N = 32 systems is rather small at the higher density, and it cuts through the second peak of the oxygen–oxygen radial distribution function, which means that the description of the second solvation shell around a given water is incomplete for the higher density system for N = 32. This inappropriate description of the second solvation shell likely affects the overall structure and dynamics of the system. In fact, the issue of the system size could be even more important for supercritical water, where large density fluctuations would require larger boxes. Thus, it would be worthwhile to extend the current study to larger systems, possibly with N = 128, 256, and beyond, to verify the convergence of the results with respect to the system size and also to examine how does the convergence of supercritical water compare with that of ambient water with respect to the system size. In the current paper, we have still chosen to present the results of N = 32 and 64 systems because these results show the extent of system size effects in this smaller region of N, which would be useful for future studies. Also, for a given N (=32 or 64), the current simulations provide us an important information that the dispersion corrections in the structure and dynamics of supercritical water are not as important as they were found for ambient water. Regarding the second point of the use of only BLYP-D functional which falls under the scheme of the so-called DFT-D2,[34,35] we note that there are many other alternative and possibly more accurate schemes available to include dispersion corrections such as DFT-D3,[64] DFT-vdW,[65] and so on. While DFT-D2 has proven to be accurate enough in many situations, the case of supercritical water can be challenging in terms of quantitative accuracy because of large density fluctuations. Certainly, it would be interesting and also worthwhile to consider some of the other schemes of dispersion corrections in ab initio simulations of supercritical water. We hope to address some of these issues in our future work.

Computational Details

The ab initio simulations were performed by using the method of Car–Parrinello molecular dynamics,[66] where the quantum many-body potentials and forces are calculated on-the-fly through quantum electronic structure calculations. The electronic structure of the extended simulation systems was represented by the Kohn–Sham (KS)[67] formulation of density functional theory within plane wave basis. The core electrons were treated via the norm-conserving pseudo potentials of Troullier–Martins,[68] and the plane wave expansion of the KS orbitals was truncated at a kinetic energy cut-off of 70 Ry. We have used a fictitious orbital mass of μ = 800 a.u.[69] and a simulation time step of 3 a.u. The hydrogen atoms were assigned the mass of deuterium to decouple the dynamics of the ionic and electronic subsystems over the simulation run length. Thus, we have considered D2O rather than H2O in the current simulations. The BLYP-D functional has been used, which incorporates dispersion corrections at the D2 level[34,35] into the BLYP[32,33] functional. As described below, we have performed simulations at fixed densities in NVT and NVE ensembles. Thus, adjustment or equilibration of the density due to dispersion corrections was not considered in the present study. It may, however, be noted that recent simulations of ambient water in the NPT ensemble have shown that the BLYP functional with D2 dispersion correction predicts the density of liquid water in close agreement with the experimental value at normal temperature and pressure.[70]

We have run a total of four different simulations with the BLYP-D functional and another four with the BLYP functional. We used the CPMD code[71] for all the present simulations. The first three systems are at density of 1.1, 0.78, and 0.39 g cm–3. All these systems contain 32 water molecules periodically replicated in all dimensions. The length of the box is 9.8652 Å for the density of 1.1 g cm–3, 11.12 Å for 0.78 g cm–3, and 13.98 Å for the density of 0.39 g cm–3, respectively. Each system was equilibrated for approximately 10 ps in a canonical (NVT) ensemble. Subsequently, each system was further run for ∼30 ps in the microcanonical (NVE) ensemble for calculations of various structural and dynamical properties. Also, to see the effects of system size, we have also run simulations of a bigger system of 64 water molecules using the BLYP and BLYP-D functionals. The edge length of the cubic box of this bigger system is 12.418 Å. These systems were equilibrated for ∼15 ps in the NVT ensemble, and then, they were run for another 30 ps in the NVE ensemble for calculations of the structural and dynamical properties. We have defined our simulation systems as systems 1, 2, and 3 for the densities of 1.1, 0.78, and 0.39 g cm–3 and N = 32, where N is the number of water molecules. System 4 is for the density of 1.1 g cm–3 and N = 64. The functional is mentioned explicitly in each case whenever the results are discussed. The BLYP results for the smaller systems of the current study are found to be in agreement with the earlier results of ref (29). The BLYP-D results at room temperature are obtained from the trajectory of ref (41) for bulk water molecules. We have verified through an additional simulation that the results of bulk water of ref (41) are very close to that of pure water for the same system size and temperature.

After the trajectories were generated, we calculated the fluctuating frequencies of OD bonds by using the wavelet method.[72−75] More details on the OD frequency calculations using this method are available in refs;[29,54,76] hence, here we include only the key details. A wavelet can be visualized as a small ripple or a wave with a compact support with an associated frequency. The wavelet transformation of a time-dependent function f(t) is defined aswhere is the complex conjugate of a mother wavelet with two tunable parameters a and b. Clearly, the right hand side of the above equation represents the overlap between the mother wavelet and the function f(t). Following the earlier work,[76] we have taken f(t) as the time-dependent fluctuations in the OD bond distance and the corresponding momentum along the simulation trajectory. We have considered the commonly used Morlet–Grossman wavelet[72] as the mother waveletwith λ = 1 and σ = 2. The wavelet transformation of the time-dependent fluctuations is then carried out for different values of a for different given values of b. The value of a, which maximizes the wavelet transform of the above equation, gives the corresponding frequency (obtained from the parameter a) at the time given by the parameter b.