Accessing the Accuracy of Density Functional Theory through Structure and Dynamics of the Water-Air Interface.

Density functional theory-based molecular dynamics simulations are increasingly being used for simulating aqueous interfaces. Nonetheless, the choice of the appropriate density functional, critically affecting the outcome of the simulation, has remained arbitrary. Here, we assess the performance of various exchange-correlation (XC) functionals, based on the metrics relevant to sum-frequency generation spectroscopy. The structure and dynamics of water at the water-air interface are governed by heterogeneous intermolecular interactions, thereby providing a critical benchmark for XC functionals. We find that the XC functionals constrained by exact functional conditions (revPBE and revPBE0) with the dispersion correction show excellent performance. The poor performance of the empirically optimized density functional (M06-L) indicates the importance of satisfying the exact functional condition. Understanding the performance of different XC functionals can aid in resolving the controversial interpretation of the interfacial water structure and direct the design of novel, improved XC functionals better suited to describing the heterogeneous interactions in condensed phases.

Density functional theory-based molecular dynamics (DFT-MD) simulation is a technique for monitoring molecular motions based on the forces calculated using DFT exchange–correlation (XC) methods.[1] While the majority of DFT-MD simulations have employed the generalized gradient approximation (GGA),[2,3] the recent increase in computational resources and efficient algorithms have allowed for the use of higher-level DFT functionals such as meta-GGA[4,5] and hybrid-GGA[6,7] in DFT-MD simulations. Because DFT-MD can deal with complex liquid–solid and liquid–gas interfaces without any empirical force field modeling, it is increasingly being used for gaining molecular-level insight into the structure and dynamics of water at aqueous interfaces.[8−14]

Within the DFT framework, the accuracy of the predicted water properties is sensitive to the adopted XC approximations. Properties of bulk water as well as interfacial water predicted by DFT by applying low levels XC approximations such as GGA are sometimes unphysical. For example, GGA functionals tend to underestimate both the density of bulk water[15,16] and its surface tension.[17,18] These limitations arise from three significant drawbacks: (1) GGA methods fail to recover the nonlocal correlation necessary to account for van der Waals (vdW) interactions.[16,19,20] (2) GGA methods, which depend only on electronic density n(r) and its gradient ∇n(r), provide inaccurate energetics, in particular, for strongly correlated systems.[21] (3) GGA methods suffer from self-interaction error.[22] These three drawbacks may be circumvented by the addition of van der Waals corrections, the extension of GGA to meta-GGA by adding a ∇2n(r)-dependency in the XC functional, and the extension of GGA to hybrid-GGA by mixing exact-exchange with GGA exchange energy, respectively. Although the accuracy of the GGA, meta-GGA, and hybrid-GGA functionals for water has been assessed in the gas phase or in the bulk, there is no rigorous and systematic assessment of the XC functionals for describing the structure and dynamics of interfacial water. Because a water molecule experiences heterogeneous intermolecular interaction at the interface unlike in the bulk, the systematic comparison of DFT-MD data at the different levels of theory can provide a unique and critical platform for examining DFT accuracy.

In this work, we explore the effect of meta-GGA and hybrid-GGA functionals on the structure and dynamics of the interfacial water together with the bulk water. The metrics for interfacial water are relevant to the sum-frequency generation (SFG) spectroscopy, allowing for direct comparison of the simulated data with the experimental results. We find that the revPBE0-D3(0) hybrid-GGA functional shows the best performance among the testified DFT methods, while the M06-L-D3(0) meta-GGA functional shows unexpected poor performance. Subsequently, we discuss the quality of the simulated SFG spectra of O–D stretch mode in isotopically diluted water, by linking with the ranking of the XC functionals, and provide insights into the controversial interpretation of the SFG spectra.[8,23−26]

To evaluate the different XC functionals, we consider the metrics for the interfacial water including the thickness of the interfacial region (δ), the fraction of the interfacial water molecules with a free O–D group, the lifetime of a free O–D group (τs), average angle of the free O–D group and the surface normal (⟨θ⟩), together with the bulk water density (ρ) and radial distribution function (RDF) of bulk water (see the details in the Supporting Information). The quantities related with the interface can be connected with various SFG measurements: The fraction of interfacial water with free O–H group was estimated to be 20–25% from the SFG measurements.[27,28] Time-resolved SFG measurement identified the lifetime of a free O–H group as 1.1 ps.[29] Furthermore, the information on the free O–H angle was obtained from the polarization-dependent SFG measurement, providing ⟨θ⟩ = 63°.[30] Although the data can be essentially obtained for the free O–D group as well, it is currently not available. Thus, we used the POLI2VS data for the reference, because the POLI2VS reproduces the data for the free O–H group accurately.[30,31] Note that the simulations did not include the nuclear quantum effects, whose effects were examined with force field-based classical MD and quantum mechanical partially adiabatic centroid MD (PA-CMD) simulations (see the Supporting Information).[32,33] Here, characterizing the free O–D group of the interfacial water is the key for computing the quantities related to the interfacial water. We used the optimal free O–D group definition.[31] More details of the calculations are provided in the Supporting Information.

The data for selected XC functionals is summarized in Table . First, we focus on bulk data. The values of the bulk H2O density in our DFT-MD simulations for M06-L-D3(0), SCAN, B97M-rV, and revPBE0-D3(0) are in good agreement with previous reports,[21,34−36] while 0.98 g/cm3 at both the B3LYP-D3(0) and HSE06-D3(0) have not been reported. Compared with the density of water predicted by the GGA functionals, the density predicted by hybrid-GGA functionals is relatively reduced, yet the density predicted by meta-GGA functionals is relatively elevated. Overestimation of the density of water with the meta-GGA functionals can be linked with the RDF data. As evident from the absence of the first minimum (hmin) in the RDF for B97M-rV and M06-L-D3(0) functionals, the hydration structure simulated with the meta-GGA functionals tends to be more disordered than that predicted by the hybrid-GGA. The improved density by the SCAN functional is attributed to the capture of dispersion forces beyond short-range. Under its influence, the water molecules at intermediate range on the H-bond network are attracted to each other by the nondirectional vdW interactions. As a result, the predicted water structure is softened by the increased population of interstitial water molecules.[35]

Table 1

Bulk and Interfacial Water Data Using Various DFT Methodsa

The average error bars for ρ, hmax, rmax, hmin, rmin, δ, Free O–D %, ⟨θ⟩, and τs, are 0.01 g/cm3, 0.05, 0.01 Å, 0.03, 0.03 Å, 1.0%, 1.0°, and 0.1 ps, respectively.

hmax (rmax) and hmin (rmin) refer to the height (position) of the first maximum and first minimum of the RDF, respectively.

Simulations performed using 128 D2O molecules.

RDF data from ref (37).

O–H data from refs (27 and 28).

O–H data from ref (30).

O–H data from ref (29).

Subsequently, we focus on the structure of the interfacial water. Our data on the thickness of the interfacial water region (δ) indicates that the meta-GGA functionals predict smaller thickness than the GGA functionals, while the hybrid-GGA functionals predict larger thickness. The opposing predictions of the meta-GGA and hybrid-GGA functionals on the thickness δ are consistent with the predicted properties of the bulk water. The different trend is also reflected in the predicted fraction of the interfacial water molecules with free O–D groups; the fraction predicted using hybrid-GGA functionals is greater than that predicted using the GGA functionals, leading to an improved description by reducing the average deviation relative to the reference from 7% to 3%. On the other hand, the meta-GGA functionals predict free O–D angles very similar to those predicted using the GGA functionals. The free O–D angle also shows the different impact of the meta-GGA and hybrid-GGA. The hybrid-GGA functionals show a deviation of 5° from the reference value, similar to the deviation for the GGA functionals, while for the meta-GGA functionals the deviation is reduced significantly to 3°. As such, the meta-GGA and hybrid-GGA functionals have a different impact on the predicted structure of water at the water–air interface.

Now, we focus on the dynamics of interfacial water. Our data show that the GGA and hybrid-GGA functionals predict the free O–D lifetimes of 1.06–1.87 ps and 0.99–1.69 ps, respectively, while the meta-GGA XC functionals predict 0.57–0.86 ps. Compared with the reference value of ∼1 ps, the meta-GGA functionals predict quite fast free O–D dynamics, while the GGA and hybrid-GGA predict relatively slower dynamics. These different time scales of the dynamics further illustrate the different impact of the meta-GGA and hybrid-GGA functionals on the predicted behavior of the interfacial water and demonstrate that a smart combination of the meta-GGA and hybrid-GGA may substantially improve the description of the interfacial water.

By using the above-listed data, we rank the performance of various XC functionals. The score for the error is computed viawhere κ is the score of the DFT method j for the calculated property i. χ denotes the value of the quantity i computed with the method j; χref is the reference value of the quantity i, and σ is the standard deviation for an extensive set of DFT-MD trajectories.[18,38] The smaller (larger) κ means that the prediction of the XC functional is better (worse). The details of the computation can be found in the Supporting Information. Figure a shows κ, the averaged κ value over i, for each of the different functionals j. This graph demonstrates that the revPBE-D3(0), SCAN, and revPBE0-D3(0) are the best XC functionals at the GGA, meta-GGA, and hybrid-GGA levels of theory, respectively. Among these, the calculation using the revPBE0-D3(0) hybrid-GGA functional shows the best performance, but at a substantially elevated computational cost (see Figure b). The revPBE-D3(0) GGA functional provides a reasonable description of the interfacial water, at a reduced computational cost. In contrast, the M06-L-D3(0) meta-GGA functional shows rather poor performance, which is somewhat surprising, given its excellent prediction of gas-phase energetics. It is however consistent with the claim of Medvedev et al.[39,40] that the XC functionals like M06-L, which does not completely obey exact functional constraints, may produce inaccurate electronic densities. Our data indicate that an accurate description of the electronic density is more critical for bulk and interfacial water than for the gas phase.

Figure 1

(a) Direct comparison of the ability of different functionals to accurately predict water properties. The smaller (larger) score κ corresponds to better (worse) predictive power of the functional. (b) Computational cost for different DFT-XC functionals. The data are normalized by the cost of the revPBE-D3(0) GGA functional.

To connect the predicted structure of water with previously reported experimental SFG spectra,[41] we computed the SFG response of the O–D stretch mode of interfacial HOD molecules in isotopically diluted water (O–D in H2O). The calculated spectra are displayed in Figure . All the simulated spectra show a sharp positive peak at 2550–2750 cm–1 and a broad negative peak centered at 2300–2500 cm–1. A positive (negative) peak corresponds to the free (hydrogen-bonded) O–D group of the interfacial water.[28,42] Compared with POLI2VS data[43] and experimental data obtained by the Tahara group,[41] we conclude that revPBE0-D3(0) is the best for reproducing the SFG features at the isotopically diluted water–air interface. The excellent reproducibility of the vibrational spectra with the revPBE0-D3(0) functionals can also be found in the infrared spectra of the bulk water.[44] Furthermore, the noticeable differences of the spectra exists with different XC functionals. For example, the spectra calculated at the PBE-D3m(BJ), BLYP-D2, and PBE-rVV10 functionals show a small, but non-negligible positive band below 2100 cm–1. This is in line with a claim of ref (24). However, Figure clearly illustrates that all the XC functionals showing a positive 2100 cm–1 band do not reproduce the water properties accurately, implying that the presence of the positive 2100 cm–1 band may arise from the poor description of the interfacial water.

Figure 2

Simulated SFG spectra of the O–D stretch mode in H2O for various XC functionals. The POLI2VS spectrum and experimental spectrum were obtained from refs (43) and (41), respectively. Because the experimental data was obtained for O–H in D2O, the frequency of the experimental data was scaled down by 0.735[46] to convert the O–H stretch frequency to the O–D frequency. Note that a positive band below 2400 cm–1 (broken line region) in the experimental data has later been attributed to an experimental artifact of the measurement and should be absent.[25,26] Each spectrum is offset by increments of 1 for clarity. The free O–D peak top of each spectrum was normalized to 2/3. The highlights of the low-frequency regions are displayed in the three panels in the left column.

We now discuss the SFG spectra simulated with the meta-GGA functionals. These show a variety of lineshapes, but all deviate somehow from the experimental data. When we focus on the data with the SCAN functional, one can notice that the negative feature is remarkably broad. In fact, the full width at half-maximum of the negative peak is 238 ± 32 cm–1, twice larger than the experimental data of 130 cm–1.[41] Even compared with BLYP-D3(0) data of 189 ± 13 cm–1, the SCAN functional predicts a very broad feature elongated to the low-frequency region. These observations may explain the controversial assignment of the SFG spectra of H2O at the water–TiO2 interface,[8,23] where BLYP-D3(0) does not show any positive SFG peaks in the low-frequency O–H stretch region, while the SCAN functional does indicate a positive peak in this frequency region. According to the current data, remodeling these SFG spectra with accurate XC functionals such as revPBE0-D3(0) or meta-hybrid GGA functionals such as SCAN0[45] is essential for resolving the controversy.

In conclusion, we have tested the quality of various GGA, meta-GGA, and hybrid-GGA functionals for a description of the structure and dynamics of water at the air–water interface. We established the significantly distinct impacts of the extension from GGA to meta-GGA and that to hybrid-GGA on the interfacial structure and dynamics. In particular, meta-GGA functionals tend to predict faster dynamics, while hybrid-GGA functionals tend to predict slower dynamics. This indicates that an appropriate combination of meta-GGA and hybrid-GGA may improve the description of interfacial water. Among the XC functionals considered here, revPBE0-D3(0) provides the best performance, while the unconstrained M06-L-D3(0) meta-GGA functional shows poor performance. Linking our results with the poor electron density prediction of the M06-L method indicates that the structure and dynamics of water at the water–air interface highlight the importance of the accurate electron density prediction.

By combining the information on the performance of the XC functional with the simulated SFG spectra of water, we found that the poor description of interfacial water tends to generate a positive low-frequency O–D stretch band. The SCAN functional tends to elongate the negative hydrogen-bonded O–D band to the low-frequency region excessively, which differs significantly from the experimental data. Again, the revPBE0-D3(0) predicts the SFG data most accurately. Our observations clearly guide the choice of the XC functional for simulating aqueous interfaces.

Computational Details

Born–Oppenheimer MD (BOMD) simulations were run using BLYP, PBE, and revPBE GGA functionals; M06-L and B97M-rV meta-GGA functionals; and B3LYP, revPBE0, and HSE06 hybrid-GGA functionals with the CP2K code,[47] while Car–Parrinello[48] MD (CPMD) simulation were run using the SCAN meta-GGA functional with the Quantum Espresso code.[49] For hybrid-GGA functionals, we used the auxiliary density matrix method (ADMM).[50] We also used various vdW correction schemes.[51−56] For describing the core electrons, we used Goedecker–Teter–Hutter pseudopotentials[57] and Hamann–Schlüter–Chiang–Vanderbilt norm-conserving pseudopotentials[58] in BOMD and CPMD simulations, respectively. We simulated 160 D2O molecules in the simulation box (L, L, L) = (16.63 Å, 16.63 Å, 44.10 Å) for BOMD simulations, while 128 D2O molecules were simulated in the box (L, L, L) = (12.44 Å, 12.44 Å, 50.00 Å) for CPMD simulation. The water–air interface is parallel to the xy plane, and the surface normal forms the z-axis. We used the NVT ensemble with the target temperature of 300 K. Furthermore, we performed the MD simulation with the POLI2VS force field model[59] for D2O.

For the SFG spectral calculations, we have used the surface-specific velocity–velocity correlation function formalism,[60] which was developed for an efficient calculation of the SFG spectra with a limited length of the trajectories. The SFG spectra of isotopically diluted water can be computed from the DFT-MD trajectories of D2O by neglecting the intra/intermolecular couplings.[60] The computed spectra are known to show higher frequency than the experimental data for the O–H (O–D) stretch mode, as the MD simulation with classical nuclei cannot account for nuclear quantum effects. To compensate the red-shift due to the nuclear quantum effects, we multiplied the computed frequency by a factor of 0.96.[61] A similar red-shift was also confirmed through the comparison of classical and quantum simulation of the SFG spectra.[62]