Vibrational mode frequency correction of liquid water in density functional theory molecular dynamics simulations with van der Waals correction

The frequencies and spectral lineshapes of the stretch and bending modes of water provide invaluable information on the microscopic structures of water in aqueous solutions and at the water/solid interfaces. Density functional theory molecular dynamics (DFT-MD) simulation has been used not only for predicting the properties of water but also for interpreting the vibrational spectra of water. Since the accuracy of the DFT-MD simulations relies on the choice of the exchange–correlation functionals and dispersion correction schemes employed, the predicted vibrational spectra at different levels of DFT theory differ significantly, prohibiting precise comparison of simulated spectra with experimental data. Here, we simulate the vibrational density of states for liquid heavy water based on various DFT-MD trajectories. We find that DFT-MD simulations tend to predict excessive inhomogeneous broadening for the stretch mode of water. Furthermore, we develop a frequency correction scheme for the stretch and bending modes of liquid water, which substantially improves the prediction of the vibrational spectra.


Introduction
The molecular structure of water is simple, yet the ensemble of water molecules leads to unique properties of liquid water, such as the higher density in the liquid phase than in the solid phase and large surface tension of water. These unique properties arise from the complex hydrogen-bond (H-bond) network of water. Understanding the relationship between the local conformation of water and H-bond strength is essential for clarifying the properties of various aqueous systems.
Information on the H-bond of water has been gained through vibrational spectroscopy such as infrared (IR) and Raman spectroscopy. 1,2 Typically, the IR spectrum of liquid water consists of three major bands; the O-H stretch mode, H-O-H bending mode, and librational modes of water. The O-H stretch and H-O-H bending modes have been widely used as reporters of the H-bond strength. [3][4][5][6] When the H-bond of water becomes stronger, the frequency of the O-H stretch mode is lowered 4 and the frequency of the H-O-H bending mode is elevated. 7 As such, frequency shifts of the O-H stretch mode in aqueous solutions reflect how the solute-water interactions alter the strength and structure of the H-bond network of water. [7][8][9][10] The target of the vibrational spectroscopy of water can be further expanded, when the vibrational technique is combined with a surface-specific technique. [11][12][13] As such, vibrational spectroscopy becomes a valuable technique for probing the microscopic structure of interfacial water.
There are, however, several challenges associated with interpreting the vibrational spectra of water. First of all, the vibrational bands of the water spectra are often broad and featureless, making it difficult to assign specific features to the molecular conformations. Secondly, for more complex, composite aqueous solutions, it is often not clear which interaction is leading to which feature of the vibrational spectra. Simulating vibrational spectra based on the trajectories sampled by molecular dynamics (MD) simulation is a direct and powerful route to interpret the vibrational spectra, because it allows us to identify the contribution of water molecules in the bulk and near the solute molecules. 5,[14][15][16][17][18] When a solute-water interaction is complicated, developing accurate force field models for computing the vibrational spectra is not straightforward. In contrast with MD simulations a Hefei National Laboratory for Physical Sciences at the Microscale, Collaborative using force field models, a density functional theory (DFT)-MD simulation computes the forces acting on the atoms based on DFT, which allows us to simulate complex solute-water interactions without using force field models. [19][20][21][22][23][24][25][26][27][28] However, the predictive powers of DFT-MD simulations regarding the water properties rely heavily on the choice of the exchange-correlation (XC) functionals and van der Waals (vdW) corrections. [29][30][31][32][33][34][35] Likewise, the peak frequencies in the vibrational spectra of liquid water are strongly affected by the XC functionals and vdW corrections. [36][37][38] To correct the peak frequency, one can use the dataset of scaling factors that are available for many XC functionals. [39][40][41] These datasets, however, do not work well for water vibrational spectroscopy, as these datasets are often created to correct for gas-phase frequencies. Furthermore, the scaling factors are generally parametrized for a set of small molecules, which may provide inaccurate numerical values for the liquid water peak frequencies. As such, generating the correction scheme for the vibrational frequency of water is urgently required for achieving reliable interpretation of the vibrational spectra of water, based on the DFT-MD simulation with available XC functionals and vdW corrections. Note that, although the vibrational spectra will be affected by several factors including nuclear quantum effects, 42,43 the scheme for correcting the frequency being presented here include these effects only implicitly; the correction factor provides a practical approach to estimate the frequency of liquid water, while accurate highlevel simulation requires nuclear quantum simulation, highlevel XC functional, and accurate basis sets, resulting in a computationally expensive simulation.
In this work, we calculated the vibrational density of states (VDOS) spectra for bulk heavy water (D 2 O) with various DFT methods and vdW methods. We find that the O-D stretch frequency is highly sensitive to the XC functionals and vdW correction schemes. In particular, DFT-MD simulations tend to predict a broader peak than the experimental data, indicating that DFT-MD describes an excessive inhomogeneous environment of bulk water. On the other hand, the D-O-D bending frequency is rather insensitive to the choice of DFT methods. Furthermore, we developed a frequency correction scheme to reproduce the stretch and bending peaks for water. These correction factors are employed to evaluate the heterogeneity of the O-D stretch mode of heavy water.
We used the mixed Gaussian and plane wave approach as implemented in the CP2K code. For the Gaussian part, the TZV2P basis set which is constructed using triple-zeta valence Gaussian basis with two sets of polarization functions was used for all BOMD simulations. We set the plane wave density cutoff of 320 Ry. Only for the M06-L-D3(0), we used a cutoff of 1200 Ry, because the M06-L functional requires a finer integration grid. 43 Norm-conserving Goedecker-Teter-Hutter pseudopotentials 67,68 were used to describe the core electrons. To accelerate the MD simulations, D 2 O was used instead of H 2 O, and the time step was set to 0.5 fs. We simulated 160 D 2 O molecules in the 16.63 Å Â 16.63 Å Â 44.10 Å cell in a slab configuration at 300 K in the NVT ensemble. The canonical sampling through velocity rescaling method 69 was employed as the thermostat. The cutoff radius of the vdW interactions was set to 10 Å.
We ran 10 independent simulations from previously generated configurations from the simulation at the revPBE-D3(0) level of theory. 32  For the SCAN meta-GGA XC functional, 70,71 we employed the CPMD methodology using the Quantum Espresso code. 72,73 We simulated 128D 2 O molecules in a 12.44 Å Â 12.44 Å Â 50 Å cell. We employed the Hamann-Schlüter-Chiang-Vanderbilt pseudopotentials 74,75 generated for PBE with a plane-wave cutoff of 85 Ry and a time step of 2 a.u. (0.0484 fs). The simulation was performed at 300 K in NVT ensemble with the Nosé-Hoover chain thermostat. 76,77 The fictitious mass of electrons was set at 100 a.u. 78,79 The 56 ps CPMD trajectory was generated after 5 ps equilibration run.
2.1.C. Classical MD. The POLI2VS model of water is known to reproduce the IR and Raman spectra of bulk liquid water as well as the SFG spectra at the water-air interface (see ESI †). [80][81][82] Based on this excellent performance of the POLI2VS model, we used the POLI2VS model of water as the reference and calculated the VDOS spectra based on the POLI2VS MD trajectories. We performed the classical MD simulation with POLI2VS model for 160 D 2 O molecules in the 16.63 Å Â 16.63 Å Â 44.10 Å cell in a slab configuration at 300 K by using the Nosé-Hoover chain thermostat 76,77 in NVT ensemble. The time step for integrating the equations of motion was 0.4 fs. The charge-charge, charge-dipole and dipole-dipole interactions were evaluated using Ewald summation, whereas quadrupole interactions were curtailed at 8.3 Å. We first ran 700 ps MD simulation for equilibration and then we obtained 1 ns for production run which was used for analysis.

VDOS spectra calculation based on MD trajectories
VDOS is computed for the bulk D 2 O molecules by using the velocity-velocity autocorrelation: where v i (t) denotes the velocity vector of the atom i at time t, m i is the mass of the atom i, T is the length of the time correlation function. We set T to 1 ps. The computed spectrum with the POLI2VS model shows response at frequencies higher than those observed in the experimental data. The higher frequencies of the POLI2VS model arise from the lack of nuclear quantum effects. 82 Thus, we scaled the frequency axis with a factor of 0.96 for both the O-D stretch and the D-O-D bending modes (see ESI †). 83 We selected the bulk D 2 O molecules as a D 2 O molecule with the z-coordinate of its oxygen atom |z| o 3 Å, where we set the origin point of z = 0 as a center of mass of the system. The region which we used for the VDOS calculation is highlighted in Fig. 1(a), while the snapshot of the simulation is drawn in Fig. 1(b).    (2) and (4), respectively. The reference data was obtained from the POLI2VS model. that the spectral shapes of the O-D stretch mode are very sensitive to the choice of the DFT functional, whereas the D-O-D bending mode is comparatively insensitive. When we focus on the linewidth of the O-D stretch peak, one can recognize that the VDOS spectra computed from the DFT-MD trajectories are generally broader than the reference VDOS spectra. For example, the VDOS spectra computed from the revPBE0-D3(0) and SCAN functionals show much broader FWHMs (299 and 327 cm À1 , respectively) than the reference POLI2VS data (FWHM of 259 cm À1 ). Since the linewidths of the spectra arise mainly from inhomogeneous broadening, this excessive broadening of the O-D stretch peak indicates that the DFT-MD simulations tend to sample a larger variety of the H-bond strengths than classical MD simulation.

VDOS Spectra
The comparison between the FWHM of the VDOS obtained from DFT-MD simulations and classical MD simulation is sharply contrasted with the comparison of their respective simulated IR spectra. Specifically, the reported linewidths of the IR stretch mode spectra using the revPBE0-D3(0) and SCAN functionals (FWHMs of 251 cm À1 36 and 246 cm À1 , 37 respectively) exhibit good agreement with the simulated IR spectra (FWHM of 257 cm À1 82 ). These different trends in the VDOS spectra and IR spectra can be explained as follows. The linewidth of the IR spectrum of neat D 2 O is determined not only by the inhomogeneous broadening but also by the intermolecular/ intramolecular vibrational coupling 1,5,84 as well as the frequencydependent transition dipole moment. 85 Thus, the good agreement between the FWHMs in the IR spectra and poor agreement between the FWHMs in the VDOS data indicates that the FWHM of the IR spectra is dominated by the inhomogeneous broadening in the DFT-MD simulations, while in the classical MD simulations, the FWHM is largely determined by intermolecular/intramolecular vibrational coupling as well as the frequency-dependent transition dipole moment.
Among the DFT-MD simulations considered in this study, the B97M-rV and M06-L-D3(0) levels of theory provide reasonable linewidths of the O-D stretch mode, in line with ref. 42. However, we also reported that these methods are inadequate for describing the water properties and the surface-specific vibrational spectra of water. 34 This means that there is room for improving the XC functionals and vdW corrections. Finally, we note that, although a small 1500 cm À1 peak can be seen in the experimental data, it cannot be properly captured in the VDOS data, as this peak arises from the combination of bending and librational motions, and the VDOS calculation cannot capture such combination band properly. A high-level computational scheme would be needed to reproduce this combination band of water. 86,87

Frequency correction factor
To correct the bandwidth as well as the peak position, we developed a correction scheme for the O-D stretch mode of D 2 O. For correcting the frequency, we assumed the expression as; where o str,ref is the characteristic frequency of the reference VDOS spectrum obtained from POLI2VS model and o str,DFT-MD is the VDOS spectrum frequency obtained from the DFT-MD trajectory. We used two parameters, a and b, because the O-D stretch frequency varies largely and thus a single scale factor is insufficient to capture the trend of the stretch mode of D 2 O. 88 To obtain the characteristic frequency, we computed the first moment of the VDOS; where o 1 and o 2 are the frequencies at which the VDOS intensities are 10% of the maximum intensity of the O-D stretch peak. The data is displayed in Table 1.
To determine two parameters, a and b, one needs to have the additional relation of o str,ref with o str,DFT-MD . This relation can be obtained from the peak frequency of the free O-D stretch in the DFT-MD simulated and experimental sum-frequency generation (SFG) spectra. At the free surface of water, a small fraction of water molecules finds itself with one non-hydrogen bonded OH/OD group that sticks into the vapor phase. The frequency of this free OH/OD group can be determined using sum-frequency generation spectroscopy. The free O-D stretch frequency of 2714 cm À1 is obtained from ref. 89, while the simulated free O-D stretch frequencies using DFT-MD simulations were obtained from ref. 34. The results are summarized in Table 1. Table 1 The simulated characteristic frequencies for the O-D stretch mode and peak frequencies for the D-O-D bending mode in the VDOS spectra as well as the free O-D stretch frequency in the SFG spectra 34 From the relation for the VDOS O-D stretch data and the free O-D stretch data, we obtained the parameters a and b. The parameters are listed in Table 2, while the O-D stretch spectra modified by using eqn (2) are plotted with the broken lines in Fig. 2. The frequency-corrected VDOS spectra of the O-D stretch mode show close agreement with the reference data. This illustrates that the frequency correction scheme is powerful for predicting the vibrational frequency and lineshape of the O-D stretch mode.
We further obtained the correction factor c for the D-O-D bending mode via; where the bending mode peak frequency in the VDOS spectra with the DFT-MD simulation is o bend,ref,DFT-MD and the reference POLI2VS VDOS spectrum is o bend,ref . The obtained parameters are also given in Table 1. In contrast with the O-D stretch mode, the peak frequency and the lineshape of the simulated D-O-D bending mode spectra are insensitive to the choice of specific DFT method and conform with the reference spectra. These correction factors can be also used to predict the vibrational spectra of H 2 O by scaling the vibrational frequency axis with the factor of 1/0.735 in the D 2 O spectra. 89 In fact, the vibrational spectra lineshapes of H 2 O and D 2 O become quite similar after the frequency axis is scaled by this factor. 90 Finally, we note that our correction scheme include a red-shift induced by nuclear quantum effects. For the DFT-MD simulation with the nuclear quantum effects such as the ring-polymer MD simulation 91 and centroid MD simulation, 92 we recommend using the coefficients, a and c, divided by a factor of 0.96.

Decomposition of the O-D stretch spectra
By establishing the frequency correction scheme for the DFT-MD trajectories, we can compare the individual spectra simulated at different DFT levels. Here, we investigate the level of heterogeneity of the O-D stretch spectra at different DFT levels of theory. To do so, we decomposed the spectra based on  the number of donating and accepting hydrogen bonds (DDAA, DDA, DAA, DD, DA, AA, and others). The hydrogen bond formation was judged at time t = 0 in eqn (1) via the geometry of water dimer; when the OÁ Á ÁO distance is less than 3.5 Å and the OÁ Á ÁO-D angle is less than 30 degrees, the hydrogen bond is formed. 93 A few representative frequency-corrected DFT-MD spectra with O-D stretch decomposition are shown in Fig. 3(a)-(d). The computed spectra with other DFT methods are given in the ESI. † The data show that the O-D stretch mode is dominated by the contribution of the DDAA species, while the other water species, such as DDA and DAA, also contribute to the water spectra. The DDAA species constitute the lowest frequency component, and the water species with a smaller number of the hydrogen bond acceptor or donor provide higher frequency component, as is expected. 4 However, one can see a marked difference of the frequency shift between the DDAA and the DDA/DAA/DD/DA/AA species. To quantify this, we computed the quantities of where o i,DFT is the center of mass frequency via eqn (3) for i = DDA, DAA, DD, DA, AA. We set Do i,ff and Do i,sd as the frequency differences, and the standard deviation for the POLI2VS data computed from DFT-MD data, respectively. The average of k i,DFT is shown in the pie chart of Fig. 3(e). This clearly indicates that one can see the marked variation of the heterogeneity of O-D stretch spectra. The BLYP and PBE levels of theory tend to predict homogeneous water than the classical MD simulations. Among these, the description highly depends on the functionals. The HSE06-D3(0), B3LYP-D3(0), revPBE-D2, revPBE-D3(BJ) and BLYP-D2 methods show small deviations from the reference POLI2VS data.

Conclusion
We have computed the VDOS spectra of heavy water for the DFT-MD trajectories with various DFT methods and compared them with the reference data. We found that different DFT methods provide a variety of the O-D stretch mode peak frequency and linewidth in the VDOS spectra. In particular, the VDOS spectra of the DFT-MD simulations tend to show a broader linewidth than the reference POLI2VS VDOS spectrum, indicating that the DFT-MD simulations excessively sample the inhomogeneous environments of water stretch mode. Subsequently, we developed the frequency correction scheme for the stretch mode and bending mode. Our data shows that the correction proposed here can improve the prediction and interpretation of the vibrational spectra of D 2 O. By further scaling the frequency axis of the O-D stretch mode to the O-H stretch mode via a factor of 1/0.735, one can also predict the O-H stretch spectra reasonably. The correction scheme was further applied to evaluate the homogeneity of the O-D stretch spectra. By decomposing the corrected VDOS into the contributions from differently hydrogen-bonded water species, we revealed that the HSE06-D3(0), B3LYP-D3(0), revPBE-D2, revPBE-D3(BJ) and BLYP-D2 levels of theory predict a reasonable homogeneity of water frequencies with the classical MD simulation.

Conflicts of interest
There are no conflicts to declare.