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.


* S Supporting Information
ABSTRACT: 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. D ensity 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 selfinteraction 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 ∇ 2 n(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  1. First, we focus on bulk data. The values of the bulk H 2 O 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/cm 3 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 (h min ) 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 shortrange. 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 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     (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 lowfrequency regions are displayed in the three panels in the left column.

The Journal of Physical Chemistry Letters
Letter By using the above-listed data, we rank the performance of various XC functionals. The score for the error is computed via where κ i j is the score of the DFT method j for the calculated property i. χ i j denotes the value of the quantity i computed with the method j; χ i ref is the reference value of the quantity i, and σ i is the standard deviation for an extensive set of DFT-MD trajectories. 18,38 The smaller (larger) κ i j means that the prediction of the XC functional is better (worse). The details of the computation can be found in the Supporting Information. Figure 1a shows κ, the averaged κ i j 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 1b). 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.
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 H 2 O). The calculated spectra are displayed in Figure 2. 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 1 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.
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 H 2 O at the water−TiO 2 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 hydrogenbonded 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−Schluẗer−Chiang−Vanderbilt norm-conserving pseudopotentials 58 in BOMD and CPMD simulations, respectively. We simulated 160 D 2 O molecules in the simulation box (L x , L y , L z ) = (16.63 Å, 16.63 Å, 44.10 Å) for BOMD simulations, while 128 D 2 O molecules were simulated in the box (L x , L y , L z ) = (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 D 2 O.
For the SFG spectral calculations, we have used the surfacespecific velocity−velocity correlation function formalism, 60 which was developed for an efficient calculation of the SFG The Journal of Physical Chemistry Letters Letter spectra with a limited length of the trajectories. The SFG spectra of isotopically diluted water can be computed from the DFT-MD trajectories of D 2 O 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 ■ ASSOCIATED CONTENT

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.9b01983.
Details of the MD and partially adiabatic centroid MD (PA-CMD) simulations, calculation protocols for target quantities, data set for MD simulations, ranking procedure, effect of ADMM, relative computational cost, simulation protocols for SFG spectra, and FWHM for the negative SFG feature (PDF)