Ab initio molecular dynamics study of collective excitations in liquid H$_2$O and D$_2$O: Effect of dispersion corrections

The collective dynamics in liquid water is an active research topic experimentally, theoretically and via simulations. Here, ab initio molecular dynamics simulations are reported in heavy and ordinary water at temperature 323.15 K, or 50$^\circ$C. The simulations in heavy water were performed both with and without dispersion corrections. We found that the dispersion correction (DFT-D3) changes the relaxation of density-density time correlation functions from a slow, typical of a supercooled state, to exponential decay behaviour of regular liquids. This implies an essential reduction of the melting point of ice in simulations with DFT-D3. Analysis of longitudinal (L) and transverse (T) current spectral functions allowed us to estimate the dispersions of acoustic and optic collective excitations and to observe the L-T mixing effect. The dispersion correction shifts the L and T optic (O) modes to lower frequencies and provides by almost thirty per cent smaller gap between the longest-wavelength LO and TO excitations, which can be a consequence of a larger effective high-frequency dielectric permittivity in simulations with dispersion corrections. Simulation in ordinary water with the dispersion correction results in frequencies of optic excitations higher than in D$_2$O, and in a long-wavelength LO-TO gap of 24 ps$^{-1}$ (127 cm$^{-1}$).


Introduction
Collective dynamics in simple and molecular liquids is so far well understood only on macroscopic length and time scales. Experiments on Brillouin scattering of light on liquids can be very nicely described by the hydrodynamic theory, which explains the main relaxation and propagation processes contributing to the measured scattered intensity of light [1,2]. However, hydrodynamic theory, which is a collection of fundamental local conservation laws (balance equations) for a system treated as continuum, is incapable of describing the atomic-scale dynamics and the collective processes on nanoscale where the atomistic structure of matter does not make it possible to treat the system as continuum. Atomistic molecular dynamics (MD) computer simulations is a very efficient tool in exploration of the dynamic properties of condensed matter and provide a lot of precious information on time-dependent correlations in liquids on nano-and atomic-scale resolution. However, the classical computer simulations make use of effective interparticle potentials, which often do not permit, especially in the case of molecular liquids, to recover the values of experimental melting points, T m , or transport coefficients. In the case of water, the difference between the calculated melting temperature from classical MD simulations and the experimental melting point of ice I h can be over 80 K (with the SPC model) [3]. Depending on the values of effective ionic charges, massless charges used for correction of a dipole moment, fixed intramolecular O-H bond length and the H-O-H angle various water models result in different melting points which are usually lower than the experimental value [3].
Ab initio molecular dynamics (AIMD) simulations, which use the density functional theory (DFT) for an electronic subsystem and classical equations of motion for ions, are free of effective interatomic potentials because the input interactions for simulations are electron-ion pseudopotentials. On-the-fly minimization of total energy to the Born-Oppenheimer surface and the consequent equilibrium AIMD runs make it possible to study instantaneous distributions of electron density that contain all the information about the bonding and polarization fluctuations in a system. The AIMD of water and ice also face some problems from the point of view of the correct melting point of ice I h . Recently, two-phase ice/water AIMD simulations with Becke-Lee-Yang-Parr (BLYP) [4,5] exchange-correlation functional indicated that at ambient pressure, the melting point of ice I h should be more than 400 K, but taking into account the socalled (London) dispersion corrections to BLYP, the melting temperature was lowered to about 360 K [6]. Namely, one of the problems of DFT calculations with generalized gradient approximations (GGA), such as BLYP, is that DFT-GGA does not reproduce the long-range attractive interaction (typical of Lennard-Jones effective potentials), whose leading term in atom pairs is ∝ −r −6 . The semiempirical dispersion correction proposed by Grimme [7] to account for van der Waals interactions in DFT-GGA calculations resulted in good agreement with experimental results in binding energies and intermolecular distances for a variety of molecular complexes. Another approach to account for van der Waals corrections in DFT calculations was initially proposed by Dion et al. [8]. Such corrections were used to study the effect of dispersion correction on the structure and self-diffusion in water [9].
One can also mention another direction in attempts to improve static and dynamic properties of water from AIMD through application of different hybrid density functionals by making use of short-ranged Hartree-Fock exchange [10], or the recent truly ab initio MD simulations using the MP2 perturbation theory for correlations beyond the Hartree-Fock energy [11]. A newer generation (D3) of the dispersion corrections was suggested in 2010 [12]. No AIMD simulations have been reported so far to check how the D3 dispersion correction would affect the value of the calculated melting point of ice I h .
Furthermore, no analysis yet regarding the time correlation functions and regarding the collective dynamics of water from AIMD studies in which dispersion corrections would be used in simulations. The collective dynamics of water itself is not completely understood, although many features such as the "fast sound" phenomenon [13,14], rotational and structural relaxation in water [15], and the gap existing between the longitudinal and transverse optic modes have been known for years. To date, there has been no generalized hydrodynamic approach to the collective dynamics in water capable of recovering the simulation results with full account for rotational and structural relaxations, although some successful attempts were made mainly to describe polarization and dielectric relaxation in water [16][17][18][19]. As concerns heavy water D 2 O, only a few studies were devoted to the analysis of collective dynamics [20][21][22]. In reference [21], the findings for heavy water were compared to the previous experimental scattering studies on H 2 O and the existence of a very flat dispersion curve at ∼ 6 meV was attributed to optic modes.
A very interesting finding for heavy water was reported from classical MD simulations using SPC/E model [22]: the analysis of longitudinal (L) and transverse (T) current spectral functions C L/T (k, ω), calculated directly from the trajectories and the corresponding velocities of particles, showed a striking feature which later on was referred to as a mixing between L and T dynamics in water: Starting from wave numbers ∼ 0.5 Å −1 , the L-and T-current spectral functions showed the emergence of a shoulder at a frequency corresponding to the T-and L-excitations, which become well defined peaks at higher wave numbers. These peaks gave evidence of both (L and T) excitations contributing to L-and T-current spectral functions. Later on, the inelastic X-ray scattering (IXS) experiments on water [23] completely supported the MD findings on the observation of the L-T coupling in reference [22]. Although there is a clear hint that rotational motion of molecules can cause the L-T coupling, so far there has not appeared either any derivation of a generalized hydrodynamic theory for water with coupled equations for L and T dynamics capable of explaining the MD or experimental observations on L-T mixing in water. We can mention here that the issue of the L-T coupling is not solely restricted to water or molecular liquids in general: The L-T coupling effects were also observed while analysing the IXS scattering experiments on liquid Hg [24] and Sn [25], as well as reported from the observation of a two-peak structure of current spectral functions obtained in ab initio simulations of liquid Sn [26] and liquid Li at high pressures [27]. No L-T coupling effects were studied so far in water by ab initio simulations. Recently, there was a report on classical simulations of supercooled water at temperature 180 K using TIP4P model [28]. It was found from a fit of several Lorentzians that for supercooled water, transverse current spectral functions contained contributions from four collective modes, while for longitudinal dynamics, only three contributions from collective excitations were estimated. Furthermore, it was shown that for the three studied systems (with density 0.939 g/cm 3 , 1.029 g/cm 3 and 1.179 g/cm 3 ) the high-frequency speed of longitudinal acoustic modes increased from 3313 m/s to 3679 m/s. Very recently, a classical simulation study of high-frequency dynamics of water using two TIP4P/2005f (with 512 molecules) and TTM3F (with 256 molecules) water models was reported [19]. The main focus was on calculations of polarization time correlation functions and estimation of dispersion of L and T optic (LO and TO, respectively) modes. The gap between the longest-wavelength LO and TO excitations obtained at the smallest accessible wave numbers in those simulations was larger than the experimental values of LO-TO gap at the corresponding temperatures. One of the shortcomings of classical MD simulations for ionic and polar liquids with non-polarizable and even polarizable effective potentials is the neglect (for non-polarizable models) or a modelling (for polarizable models) of high-frequency polarizability of electron subsystem. It was shown in reference [29], from a comparison of classical simulations using a rigid ion model for molten salt NaCl and ab initio simulations at the same thermodynamic point, that the high-frequency polarizability leads to a "softening" of the long-wavelength LO frequency by a factor ∝ ε ∞ , while the transverse TO frequency becomes slightly higher (TO "hardening") in comparison with the non-polarizable case. Here, ε ∞ is the high-frequency dielectric permittivity. These findings were supported by sequential AIMD simulations for molten salts NaI [30], RbF [31] and LiBr [32]. Therefore, it would be interesting to study the LO-TO gap for most long-wavelength excitations in water from ab initio simulations.
In this study, we aimed at calculating, from AIMD, the L-and T-current spectral functions and estimating the dispersion of collective excitations. Initially, we intended to simulate only heavy water D 2 O both with and without dispersion corrections to DFT in order to estimate the effect of dispersion corrections on the decay of density-density time correlation functions, which is very sensitive to the regular liquid/supercooled/glassy state of the studied system, as well as to observe the effect of dispersion corrections on the high-frequency speed of sound and LO, TO frequencies. Additional simulations were performed in light, or "normal" water, with dispersion corrections applied, which enabled us to compare our results on the dispersion of LO and TO excitations with the recently reported ones [19]. Another important issue was to estimate whether the L-T coupling effects are observed in ab initio simulations of water and heavy water, and to estimate in general in what exactly consists the difference in dispersion of collective excitations between heavy and regular water.
The remaining Paper is organized as follows: In the next section we supply details of our ab initio simulations, while the third section contains results for static structure of D 2 O obtained from simulations with and without dispersion corrections; the static structure in light water is practically identical to that in heavy water. Further, we will show the results for the density-density time correlation functions, L-and Tcurrent spectral functions, and dispersions of collective excitations. The last section contains conclusions of this study.

Molecular dynamics simulations
We performed ab initio simulations having a collection of 128 molecules in a periodic box with the box length of 15.7459 Å for three systems: heavy water D 2 O with D3 dispersion correction and without it, and ordinary water H 2 O with D3 dispersion correction. The details of our simulations are very similar to those used in an earlier work [33]. We performed Born-Oppenheimer molecular dynamics simulations using the density functional theory to provide the forces into the Verlet algorithm. We integrated the equations of motion with a time step of 0.5 fs in the canonical, NV T , ensemble, and used a Nosé-Hoover thermostat chain with a time constant of 2.4 ps to equilibrate the average temperature at 323.15 K. The slightly elevated temperature is motivated by the known deficiency of the approximations in the details of the DFT treatment discussed below [33][34][35]. Either the mass of hydrogen or deuterium was used in simulations of light and heavy water, respectively. All the trajectories are at least 50 ps long, and the first 10 ps were neglected as a period of equilibration.
We used the BLYP generalized gradient approximation as the exchange-correlation term in the Kohn-Sham equations of DFT. Further, the D3 approximation [12] was added to the Kohn-Sham expression of total energy in order to incorporate the London dispersion forces missing in the BLYP, or (semi-)local approximations in general.
The CP2K suite of programs [36] was used in the simulations, in particular its QuickStep [37] module. The Gaussian Plane Wave (GPW) method [38] was used with the triple-zeta valence doubly polarized (TZV2P) basis set of Gaussians and plane waves up to 400 Ry to expand the Kohn-Sham orbitals and the electron density, respectively. We also employed the NN50 smoothing algorithm [37] to reduce the errors in the evaluation of the energy and forces of the BLYP term, since tests [33] have shown that this combination delivers a good convergence in the simulations. The action of the nuclei and core electrons on the valence electrons was replaced by the Goedecker-Teter-Hutter (GTH) norm-conserving pseudopotentials [39,40].
At each configuration, we calculated Fourier-components of the following dynamic variables: partial particle densities and densities of partial transverse momentum Here, N O and N H/D are 128 and 256, respectively, in our simulations. The smallest accessible wave number from our AIMD was 0.399 Å −1 . We studied the range of wave numbers up to 4.5 Å −1 . All the wave number-dependent quantities calculated from AIMD were averaged over all possible directions of wave vectors corresponding to the same absolute value. The partial densities were used for calculations of collective Bhatia-Thornton time correlation functions and corresponding static structure factors either in "T-C" (total number-concentration) or in "M-X" (total mass-mass-concentration) representations [41]. All these partial and collective representations contain the full information about the same collective mode, although, as it was shown in [42,43], the mass-concentration currents are the most useful ones in the studies of optic collective modes.

Results and discussion
Partial pair distribution functions from simulations with and without the D3 dispersion correction are shown in figure 1. The simulations of D 2 O and H 2 O both with dispersion correction resulted in practically identical static, or average, structure. One can see in figure 1 that the D3 correction led to a strong reduction in the depth of the first minimum of partial O-O distribution compared to the regular BLYP simulations, which resulted in an overbonded water molecules; this is related to a very high melting temperature of ice I h from AIMD with BLYP exchange-correlation functional, suggested in reference [6]. The height of all the first intermolecular peaks in the partial distribution functions was reduced in BLYP+D3 simulations, and first intermolecular minima shifted to larger values. This means that there is more active exchange of molecules in the first coordination shell with molecules from the bulk, being consistent with an enhanced diffusion constant found in reference [33]. The inclusion of the D3 dispersion correction slightly reduced the values of intramolecular O-H distance and H-O-H angle, as in the right-hand frames of figure 1 one can see that the corresponding distributions were shifted towards smaller values in simulations with BLYP+D3. In simulations with classical equations of motion for ions, there is practically no difference in the results obtained for H 2 O and D 2 O, which appears in quantum treatment via The large peak at about 105 • corresponds to tetrahedral ordering of the nearest neighbors. One can see that the application of D3 dispersion correction significantly reduces the contribution from tetrahedral configurations of nearest neighbors, which is actually in line with a suggestion of the reduction of the melting point of ice in simulations with D3 dispersion correction. In the literature one can find a similar  [46] show pronounced maxima at about 100 • − 105 • and a small shoulder close to 60 • , which is in agreement with our BLYP+D3 simulations.
Partial static structure factors are of huge importance in understanding the role of partial contributions to the measured X-ray or neutron diffraction intensities. Here, we report the Bhatia-Thornton total density S tt (k) and concentration S cc (k) static structure factors. The total density structure factor shown in figure 3 has the main peak at ∼ 2 Å −1 and a well-defined shoulder in wave number range k ∼ 2.5−2.8 Å −1 , which is in agreement with experimental results [50] and ab initio calculations of X-ray spectra of liquid water [51]. The concentration structure factor S cc (k) has a well-defined maximum at k ∼ 3.1 Å, which is related to short-range order of hydrogen bonds. The concentration structure factor should have the same long-wavelength asymptote of ∼ k 2 as the charge structure factor. In figure 3 (b), we show the longwavelength asymptotes of concentration structure factors in simulations with and without D3 dispersion correction, which nicely recover the ∼ k 2 behaviour. The coefficient at the k 2 term depends on the value of the high-frequency dielectric permittivity ε ∞ , which makes it possible to estimate ε ∞ if classical simulations with non-polarizable effective interaction models were available at the same thermodynamic  the relaxation of F tt (k, t) is typically a simple exponential one as for regular liquids. point, like it was done in reference [29] in liquid NaCl. Here, figure 3 (b) implies that the system simulated with D3 dispersion correction has a larger value of (effective) ε ∞ . This fact is important in analysing the frequencies of longitudinal optic modes that experience a softening for a system with larger values of ε ∞ .
Electron structure simulations further permit to perform calculations of more general total charge structure factors S QQ (k), where Q(k, t ) are the total charge fluctuations composed of point positive ions and instantaneous distributions of electron density in AIMD. Recently, it was shown that the long-wavelength asymptote of S QQ (k) in water recovers the ∼ k 2 behaviour [52].
The total density-density time correlation functions F tt (k, t ) at three wave numbers, calculated from simulations of D 2 O with and without D3 dispersion correction, are shown in figure 4. In general, the decay of density-density time correlation functions is very sensitive to the studied thermodynamic point on the phase diagram. At low wave numbers, the hydrodynamic theory predicts leading contributions to F tt (k, t ) coming from longitudinal collective excitations (and represented by a damped oscillator) and from thermal relaxation having an exponential tail. The exponential decay of the tail of density-density time correlation inherent to an ordinary liquid state is replaced by a slow stretched-exponential decay if the liquid becomes supercooled. By approaching the glass transition, the tail of the F tt (k, t ) shows features connected with α-relaxation, and ultimately, at the glass transition, there appears a non-decaying plateau F tt (k, t → ∞) = f (k) known as the non-ergodicity factor [53]. Moreover, for supercritical fluids it was recently shown that the analysis of density-density correlations allows one to discriminate between liquid-like and gas-like types of fluid in the supercritical regime [54,55].
In figure 4, one can see that the total density time correlation functions obtained from simulations without D3 dispersion correction show slow relaxations typical of supercooled liquids, while our AIMD simulations with D3 dispersion correction give evidence of regular exponential tails of F tt (k, t ). This is a further proof that the simulated thermodynamic point is below the melting line in simulations without accounting for dispersion corrections, while BLYP+D3 simulations result in density-density correlations typical of a regular liquid state.
Longitudinal and transverse total current spectral functions C L/T tt (k, ω) obtained from simulations of and without D3 dispersion correction in D 2 O resulted practically in the same frequency. These peaks correspond to intramolecular bending normal modes, which have nearly two times lower frequencies than those of the symmetric stretch normal modes. The frequencies of L and T acoustic collective excitations can be observed as the low-frequency peaks in the corresponding total current spectral functions C L/T tt (k, ω). In figures 5 and 6 one can see the emergence of another low-frequency peak at k > 1 Å −1 in C L tt (k, ω), which practically coincides with the frequency of the transverse excitation. This is the effect of L-T mixing reported previously from classical MD simulations [22] and X-ray scattering experiments [23].  One can easily observe the optic-like modes as peaks in the mass-concentration current spectral functions C L/T xx (k, ω), shown for light water in figure 7. A nice property of C L/T xx (k, ω) spectral function is that the contributions from acoustic modes in it are suppressed and the high-frequency optic modes are well observable, in contrast to the C L/T tt (k, ω) spectral function. The L and T optic modes clearly appear at different frequencies in figure 7, that is typical for ionic melts and polar fluids, and result in a LO-TO gap in the long-wavelength region. For ionic crystals, there is a famous Lyddane-Sachs-Teller relation for the difference between long-wavelength LO and TO phonons (non-damped propagating modes of small where ε(0) is the macroscopic static dielectric permittivity of a non-conducting system. Although in nonconducting liquids there should be some corrections to the Lyddane-Sachs-Teller due to the damping of optic modes, which in its turn renormalizes the observed LO and TO frequencies, this relation gives a hint of how the high-frequency screening by the electron density reduces the LO-TO gap.
The dispersion curves of acoustic and optic collective excitations obtained from the L and T current spectral functions for D 2 O with and without D3 dispersion correction are shown in figure 8. The small size of the simulated system does not allow one to observe the transition towards hydrodynamic dispersion law with adiabatic speed of sound. The dashed straight lines indicate the linear dispersion law with the apparent high-frequency speed of sound. We obtained a reduction in the value of the apparent highfrequency speed of sound from 3640 to 3070 m/s when the D3 dispersion correction was switched on. At the lowest accessible wave number, the value of the LO-TO gap was also lower when the D3 dispersion 23604-9 correction was applied. The longitudinal LO frequency was reduced from 123 to 107 ps −1 , while the transverse TO frequency changed from 113.8 to 102.8 ps −1 , resulting in a reduction of the LO-TO gap from 16 to 11 ps −1 due to the effect of D3 dispersion correction, which transforms the collective dynamics typical of an undercooled state to the regular dynamics of liquid as well as increased the high-frequency dielectric permittivity [see figure 3 (b)]. Both simulations with and without D3 dispersion correction show the same L-T mixing as was reported earlier from classical MD simulations [22]. In both cases, there is also a similar behaviour in the TO mode, whose frequency decreases with an increasing wave number up to k ∼ 1.5 Å −1 . At larger wave numbers, it changes the dependence, showing a small increase up to a wave number k ∼ 2.5 Å −1 , where it merges with the LO branch. At higher wave numbers, LO and TO branches stay at the same frequencies. This behaviour of dispersion implies that there is a similar crossover between the intrinsically collective (k < 1.5 Å −1 ) and partial (at k > 1.5 Å −1 ) types of dynamics described in detail in reference [42].
In figure 9, we present the dispersion of longitudinal and transverse collective excitations in liquid H 2 O from simulations using BLYP+D3 approach. The value of apparent high-frequency speed of sound is 3223 m/s, which is in good agreement with the values ∼ 3100 m/s obtained from the X-ray scattering experiments on water close to the melting point [57] and ∼ 3200 m/s reported in [21]. The optic branches of collective excitations due to the smaller molecular mass are located at higher frequencies than in the case of heavy water. At the long-wavelength limit, the LO and TO frequencies are 152.2 and 128.2 ps −1 , respectively, yielding our estimate for the LO-TO gap of 24 ps −1 . The frequency of the LO long-wavelength exci- tation is slightly lower than the one obtained from classical simulations using a flexible non-polarizable model TIP4P/2005f [19]. This difference could be due to a softening of the LO mode in AIMD due to the high-frequency screening effects. Our calculated TO frequency is higher than the one from the classical simulations [19], thus both effects resulting in a smaller long-wavelength LO-TO gap from AIMD.

Conclusions
We performed ab initio molecular dynamics simulations for heavy water with and without the third generation dispersion correction of Grimme and co-workers D3 to DFT [12]. We found that the D3 dispersion correction strongly affects the static and dynamic structure of heavy water. Very sensitive to the thermodynamic state, the behaviour of density-density time correlation functions showed slow, stretched exponential relaxation, typical of supercooled liquids, in the case of simulations without D3 dispersion correction, and an exponential relaxation of regular liquids when the dispersion correction was used. This finding could be a consequence of a significant reduction of the melting point of ice in the DFT/BLYP+D3 simulations compared to BLYP-only approach.
Longitudinal and transverse current spectral functions were analyzed in order to estimate the dispersion of acoustic and optic excitations. The application of the D3 dispersion correction resulted in a significant reduction of the apparent high-frequency speed of sound, as well as in a reduction of the gap between long-wavelength LO and TO excitations. This fact indicates that the D3 dispersion correction leads to higher values of the high-frequency dielectric permittivity ε ∞ . This change is effective, or indirect, since the D3 dispersion correction does not directly involve the electronic structure but only via the changed static and dynamic structure of the water.
Another simulation with D3 dispersion correction was performed for light water at the same temperature and number density, and it resulted in an identical static structure as in D 2 O with the D3 dispersion correction. Due to a smaller molecular mass, higher frequencies of optic L and T and of intramolecular modes were obtained. We have a good agreement for the value of the high-frequency speed of sound with X-ray scattering experiments. The gap between the long-wavelength LO and TO excitations in light water was obtained as ∼ 24 ps −1 , or 127 cm −1 . The calculated frequency of the LO excitations from our AIMD is slightly lower than the one from classical simulations using flexible non-polarizable model TIP4P/2005f [19], while frequencies of transverse optic excitations in our ab initio simulations with D3 dispersion correction were much higher than the ones from classical simulations.
In all three simulations, we observed the L-T mixing effect, since the L/T current spectral functions contained an additional peak, corresponding to the T/L excitations, at relatively large wave numbers. This L-T mixing effect was observed only for acoustic excitations; the optic modes appeared at separated frequencies at wave numbers k < 2.5 Å −1 . So far there is no indication that TO modes can appear in the longitudinal spectra.
Our study gives indications of a significant reduction of the ice I h melting temperature in ab initio simulations with D3 dispersion correction applied. There is thus a quest for large AIMD simulations with D3 dispersion correction on the melting of ice I h .