Structure and Dynamics of Water at the Water–Air Interface Using First-Principles Molecular Dynamics Simulations. II. NonLocal vs Empirical van der Waals Corrections

van der Waals (vdW) correction schemes have been recognized to be essential for an accurate description of liquid water in first-principles molecular dynamics simulation. The description of the structure and dynamics of water is governed by the type of the vdW corrections. So far, two vdW correction schemes have been often used: empirical vdW corrections and nonlocal vdW corrections. In this paper, we assess the influence of the empirical vs nonlocal vdW correction schemes on the structure and dynamics of water at the water–air interface. Since the structure of water at the water–air interface is established by a delicate balance of hydrogen bond formation and breaking, the simulation at the water–air interface provides a unique platform to testify as to the heterogeneous interaction of water. We used the metrics [Ohto et al. J. Chem. Theory Comput., 2019, 15, 595−60230468702] which are directly connected with the sum-frequency generation spectroscopic measurement. We find that the overall performance of nonlocal vdW methods is either similar or worse compared to the empirical vdW methods. We also investigated the performance of the optB88-DRSLL functional, which showed slightly less accuracy than the revPBE-D3 method. We conclude that the revPBE-D3 method shows the best performance for describing the interfacial water.


INTRODUCTION
van der Waals (vdW) interactions are essential for predicting structure and dynamics of water in the bulk and at the interface. 1,2 Although the density functional theory (DFT) itself is rigorous, the approximated exchange-correlation (XC) functionals within the generalized gradient approximation (GGA) or hybrid GGA cannot compute the vdW interaction energy accurately. Thus, first-principles molecular dynamics (FPMD) simulations performed with the GGA-XC or hybrid GGA-XC functionals erroneously predict the properties of water. For example, FPMD simulations based on the GGA-XC functionals without vdW corrections underestimated the bulk density of water. The underestimated density values of water in FPMD simulation is known to be improved largely by including the vdW corrections. 3−9 Analogously, the simulated values of the surface tensions of water at the water−air interface are drastically underestimated without the vdW corrections, whereas they are improved with the vdW corrections. 10,11 Since the impact of the vdW corrections is huge in predicting the properties of water, the choice of the vdW correction needs to be examined carefully.
So far, the frequently used vdW corrections in FPMD simulations of liquid water were the empirical vdW corrections based, in particular, on Grimme's D2 12 and D3 13−16 corrections. The D2 correction employs the environment-independent parameters for computing the vdW correction energy, while the D3 correction used the environment-dependent parameters.
For both methods, the total DFT energy of a system, E tot , within the GGA 17 (1) where n(r) denotes the electron density, T s [n(r)] denotes the kinetic energy, E ext [n(r)] denotes the electron-nuclei interaction energy, E H [n(r)] denotes the electron−electron interaction energy, E X GGA [n(r)] denotes the exchange energy of the GGA functional, E C GGA [n(r)] denotes the correlation energy of the GGA functional, and E empirical-vdW (r) denotes the vdW correction energy based on the empirical vdW method. Since the empirical vdW corrections are computationally inexpensive (less than ∼20% of the computational cost of the base XC functional), they have been widely used for computing the properties of water.
Beyond the empirical vdW corrections, a computational scheme to evaluate the vdW correction energy within the DFT scheme has been proposed. This vdW correction scheme is called the nonlocal vdW (nl-vdW) correction scheme. The total DFT energy is given as: where E C GGA/LDA [n(r)] is the correlation energy calculated within the GGA or the local density approximation (LDA), and E C nl-vdW [n(r)] denotes the correlation energy computed in the nl-vdW correction scheme. Since the nl-vdW corrections are calculated within the DFT framework through the electron density, the computational cost of GGA-FPMD simulations with nl-vdW corrections is more than doubled, compared with the GGA-FPMD simulation without any vdW corrections. The frequently used nl-vdW correction schemes are DRSLL 20,21 and rVV10. 22,23 The nl-vdW corrections were combined with various GGA XC functionals based on eq 2. 9,24−28 The DRSLL nl-vdW corrections have been used for predicting the property of bulk water. 25,26 Furthermore, the optB88-DRSLL 29 was suggested as one of the best DFT methods for computing the water properties. 1,9,24,30,31 However, the water properties in the bulk and at the interface predicted with the nl-vdW correction schemes have not been rigorously compared with those with the empirical vdW corrections.
Here, we evaluate various empirical vdW and nl-vdW corrections by focusing on the structure and dynamics of water at the water−air interface. 10,32−38 At the water−air interface, the free O−H (O−D) groups of the interfacial (heavy) water stick out, characterizing the structure of the interface. 32,33,39−44 Since the behavior of the free O−H (O−D) groups is well characterized by the surface-specific sumfrequency generation (SFG) spectroscopy technique, 45−51 we shall focus on the properties of the free O−H (O−D) groups. Furthermore, we also compute several bulk water properties. Based on the computed data, we make a ranking of the performance of the vdW correction schemes, with the similar scheme to ref 11. We find that the empirical D3 vdW correction scheme is comparably good as the rVV10 nl-vdW correction and better than the DRSLL nl-vdW correction. This indicates that the empirical vdW correction schemes are adequate to improve the description of interfacial water in FPMD simulation at the GGA level of theory.

VDW CORRECTIONS
2.1. Environment-Independent Empirical vdW Corrections. A simple description of the vdW correction employs environment-independent empirical vdW parameters. A typical vdW correction in this category is given by the Grimme's D2 model, 12 which can be written as where C 6 AB is the coefficient for atoms A and B, r AB is the distance between A and B, s 6 is a scaling factor, and f damp (r AB ) is the shortrange damping function. C 6 AB is empirically determined and is independent of the local environment, making this correction scheme categorized into the environment-independent empirical vdW correction scheme. 52 2.2. Environment-Dependent Empirical vdW Corrections. Straightforward extension of the Grimme's D2 correction is to make the coefficient C AB dependent on the local environment. The vdW correction schemes using the environment-dependent C AB (CN A ,CN B ) coefficients are called the environment-dependent empirical vdW correction, 52 where CN A and CN B stand for the coordination number for atoms A and B, respectively. In this scheme, Grimme's D3 approach 13,14,16 has been used frequently. The vdW correction energy within the Grimme's D3 scheme is given by In addition to the environment-dependent C AB (CN A ,CN B ) coefficient, the D3 model differs from the D2 model in that the D3 model contains the higher order interactions term (n = 8) and the three-body interaction term (E (3) ).
The Grimme's D3 correction has mainly three versions: D3(0), 13 D3(BJ), 14 and D3m(BJ). 15 The D3(BJ) model uses a Becke-Johnson damping function 53 for short-range interactions, while the D3(0) model uses a Chai-Head-Gordon damping function. 54 The parameter set in the D3m(BJ) model was modified from the parameter set of the D3(BJ) model, by using a larger database of the conformational energy. 15 2.3. Nonlocal vdW Density Functional Correction Schemes. The nl-vdW correction utilizes the DFT framework for computing long-range vdW interactions and without the need to define the pairwise interactions terms. The nonlocal correlational energy E C nl-vdW [n(r)] can be written as E n d d n n r r r r r r r ( ) 1 2 ( ) ( , ) ( ) where the correlation kernel ϕ is a function of the density n(r) and the gradient of the density ∇n(r). The different versions of the nl-vdW correction arise from the different functional forms of ϕ. 55−57 So far, two popular nl-vdW correction schemes were presented: DRSLL 20,21 and rVV10. 22,23 For the DRSLL method presented by Dion and co-workers, 20 the nonlocal correlation kernel ϕ is directly derived using the plasmon-pole approximation and evaluated in terms of an auxiliary function q = q(n(r),|∇n(r)|) in the DRSLL. On the other hand, the VV10 method was originally developed and revised by Vydrov and van Voorhis 22 and was optimized by Sabatini and co-workers, 23 which is called rVV10. The rVV10 utilizes a predefined empirical ϕ which has separate dependences for n(r) and ∇n(r). In the rVV10 approach ϕ contains two parametersb and C: b controls the short-range damping of asymptotic ϕ, while C controls the long-range behavior of ϕ. 21,28,56 In addition to the different functional forms of ϕ, the DRSLL and rVV10 use different correlation energies; the DRSLL uses the correlation energy calculated within the LDA (E C LDA ), while the rVV10 method uses the correlation energy calculated within the GGA (E C GGA ).
We have used a triple ζ valence Gaussian basis set with two sets of polarization functions (TZV2P). Norm-conserving Goedecker-Teter-Hutter pseudopotentials 63,64 were used to describe the core electrons. We utilized the self-consistent field

Journal of Chemical Theory and Computation
Article (SCF) cutoff of 2 × 10 −7 hartrees and kinetic energy density cutoff set at 1 × 10 −10 . The plane wave density cutoff was set to 320 Ry. 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. All simulations were performed at 300 K in the NVT ensemble with the CSVR themostat. 65 We used 160 D 2 O molecules in a simulation cell (L x ,L y ,L z ) = (16.63 Å, 16.63 Å, 44.10 Å), where the water−air interface is parallel to the xy plane and the surface normal is parallel to the z-axis. The cutoff of the vdW interactions was set to 10 Å. We ran 10 independent simulations from previously generated configurations from the revPBE-D3(0) functional for revPBE+vdW simulations and the PBE-D3(0) functional for PBE+vdW and optB88-DRSLL simulations. Fifteen ps of the production run was generated after 5 ps equilibration runs. Configurations were captured at every 10th step (i.e., at 5 fs) and were used for analysis.
3.2. POLI2VS MD Simulation. We generated the reference bulk and interfacial data from the MD trajectory with the POLI2VS model. 66 The details for the simulations have been provided previously. 11 The POLI2VS polarizable water model is the most tested water model at the water−air interface, through the direct comparison between the SFG measurement and SFG simulation. 33,66−73 Furthermore, the direct comparison of the POLI2VS model and another accurate model (MB-pol) shows that the description of interfacial water given by POLI2VS and MB-pol models is within the error bars. 33 Thus, we used the POLI2VS model to generate the reference data.

TARGET QUANTITIES
4.1. Density Profile of Interfacial Water. The bulk water density ρ 0 (at zero pressure) can be calculated by fitting the density profile in the water slab system to the hyperbolic tangent function ρ(z) along the surface normal (z-axis) i where z G is the z-coordinate of the Gibbs dividing surface, and δ is the interfacial thickness parameter. For this calculation, the center of mass of the water slab was set at the origin of the z-axis, and the deuterium atoms in D 2 O were replaced by hydrogen atoms. Note that ρ 0 calculated in the slab model agrees well with the density calculated in the NPT ensemble within the error. In fact, the water density predicted in the slab model (0.989 ± 0.001 g/cm 3 ) 11 is very close to the density predicted in the NPT ensemble (0.993 ± 0.001 g/cm 3 ) 66 in the MD simulation with the POLI2VS force field model. An O−D group of the D 2 O molecule in the interfacial region is defined as free, when the intermolecular O···O distance of its oxygen atom and oxygen atom of any other water molecule is greater than 3.5 Å and the D−O···O angle is more than 50°. 32 Otherwise, an O−D group is defined as being hydrogen-bonded. The interfacial region for the slab system is defined along the zaxis as 45 The fraction of the interfacial free O−D groups is calculated as the sum of the DA and DAA fractions, which have been defined elsewhere. 11, 32 The fraction of the free O−D group estimated from the previously generated reference POLI2VS value is 28%, which is in good agreement with SFG spectroscopy measurements of 20−25%. 45,74 We can also calculate the averaged angle ⟨θ⟩, formed by the free O−D groups and the surface normal. This value can also be extracted using polarization-dependent SFG spectroscopy. 33

RESULTS
The obtained data is summarized in Table 2

Journal of Chemical Theory and Computation
Article the density of water simulated with the D3 vdW correction agrees with the previous reports. On the other hand, the calculated density value with the nl-vdW correction has not been so much reported. The bulk water density for the PBE-DRSLL (ρ = 1.12 ± 0.01 g/cm 3 ) is similar to the reported values of 1.13 g/cm 3 in ref 26 and 1.18 g/cm 3 in ref 28. However, our data for the revPBE-DRSLL differs from the data in ref 26; we obtained ρ = 0.93 ± 0.01 g/cm 3 , smaller than the value of 1.02 g/cm 3 . 26 This can be rationalized by the different techniques to compute the density of water. We used the slab model in the NVT ensemble, while ref 26 used the computed the pressure−density relation and estimated the density at 1 atm. As is mentioned above, we have already checked that the water density predicted in the slab model is very close to the density predicted in the NPT ensemble. Together with the very good agreement with our estimation of the water density in the slab model with that in the NPT simulation in the GGA+D3 level of theory, 11 we believe that the water density at the revPBE-DRSLL level of theory is 0.93 ± 0.01 g/cm 3 . We would like to note that, to the best of our knowledge, simulated density of water using the PBE-rVV10 and revPBE-rVV10 methods has not been reported yet. Now, we turn our focus on the effects of various vdW corrections on the bulk water density. Table 2 indicates that, compared with the FPMD simulation without the vdW corrections, the density values are improved with the vdW corrections, while there is no clear trend that the nl-vdW correction scheme provides better density value than empirical vdW corrections. In fact, the PBE-DRSLL level of theory provides the largest deviation from the reference value of 0.989 g/cm 3 .

Interfacial Thickness.
For the interfacial thickness δ, the tendency is consistent with the bulk water density. The thickness parameter decreases with any type of the vdW correction schemes, making the thickness parameter closer to the reference value of 1.292 Å. On the other hand, we could not see a signature that the nl-vdW corrections provide the better thickness value than the empirical vdW corrections. Again, the largest deviation from the reference value is obtained at the PBE-DRSLL level of theory.

Fraction of Interfacial Water Molecules with Free O−D Groups.
Here, we discuss the effect of the vdW corrections on the fraction of the interfacial water molecules with the free O−D groups (DA-type and DAA-type water conformation). For the PBE and revPBE cases, the fraction of DA water molecules decreases drastically with any type of the vdW correction schemes, while the fraction of DAA is rather insensitive to the vdW corrections. Interestingly, the PBE-DRSLL drastically underestimated the DA fraction (6%), while the revPBE-DRSLL overestimated the DA fraction significantly (19%). As a result, the DRSLL correction provides noticeably scattered values of the DA fraction. The rVV10 method predicts the DA fraction closer to the reference value than the DRSLL, whereas it is worse than the empirical D3 vdW corrections. Overall, we conclude that the prediction of the free O−D fraction with the nl-vdW corrections is less accurate than that with the empirical D3 vdW corrections.
5.4. Angle of the Free O−D Group of Interfacial Water Molecules. We do not see any systematic trend on the effects of the vdW correction on the free O−D angle ⟨θ⟩ of the interfacial water molecules. However, it is worth mentioning that the worst prediction of the free O−D angle with the PBE XC functional is given by the rVV10 nl-vdW correction scheme, while the worst prediction with the revPBE XC functional is given by the

Journal of Chemical Theory and Computation
Article DRSLL nl-vdW correction scheme. As a result, we conclude that the use of the nl-vdW correction scheme does not improve the description of the free O−D angle over the empirical vdW corrections.

Lifetime of Free O−D Groups.
We computed the time correlation function for the free O−D group at the water−air interface and obtained the lifetime of the free O−D group from the fit of eq 9. Figure 1 displays the time correlation function C(t) for the DFT XC functionals, and Table 2  In our previous report, 11 the dynamics of the free O−D group is strongly affected by the XC functional rather than the vdW corrections. This trend can be seen in the current data except for the DRSLL vdW correction; τ s calculated with the DRSLL correction is almost half of τ s with the D3(0) vdW correction. This is likely due to an underestimation of the hydrogen bond strength in the DRSLL method. 80,81 Indeed, this is in line with the results of hydrogen bond dynamics in ref 9, where Klein and co-workers have reported that the lifetime of the hydrogen bond is 5.0 ps at the PBE-D3(0) level of theory and 1.0 ps at the PBE-DRSLL level of theory. As such, we conclude that the interfacial water dynamics is mainly governed by the XC functional and thus the effect of the vdW correction is minor, except for the DRSLL correction. 5.6. RDF. We plot g OO (r) in Figure 2 and summarize the heights and positions of the first maxima (h max ;r max ) and minima (h min ;r min ) in Table 2. First, we compare our simulated data with the previous reports. The RDF with the empirical D2 and D3 correction schemes has been reported in many papers, 3,9,34,82 and our data is consistent with this literature. In contrast, the RDF data with the nl-vdW corrections are available only for the DRSLL correction. For the simulation at the PBE-DRSLL level of theory, it was reported that h max = 2.77 and r max = 2. 80  The data in Table 2 indicates that the RDF data is sensitive to the XC functional and is rather insensitive to the vdW corrections, except for the DRSLL correction. The DRSLL corrections drastically reduce h max . 5.7. Performance of the optB88-DRSLL Functional. The optB88-DRSLL functional was developed by Michaelides and co-workers 29 based on the original B88 exchange functional. 19 This functional is known to predict the energy of water and ice clusters well. 1 Furthermore, RDF 9,24 and dipole moment 24 obtained from the optB88-DRSLL simulation of bulk water is in good agreement with the experimental data. The reported bulk water density in the simulation at the optB88-DRSLL level of theory was 1.08 g/cm 3 at 295 K and 1 atm, 25 which is consistent with our density value of ρ = 1.07 ± 0.02 g/ cm 3 . This density value is higher than the experimental value of ρ = 1.00 g/cm 3 . As such, the performance of the optB88-DRSLL in the bulk is good but still far from being perfect.
How about the description of the optB88-DRSLL level of theory for the water−air interface? The fraction of the interfacial water molecules with the free O−D groups is 20%, which is much lower than the reference value of 28%. The angle formed by the free O−D group and surface normal is ⟨θ⟩ = 52°, which is also drastically underestimated. Figure 2 indicates that the optB88-DRSLL predicts slower free O−D group dynamics (τ s = 1.69 ps). As such, we conclude that, within the metric defined in this study for air−water interface, the performance of the optB88-DRSLL level of theory is not so good.
5.8. Performance Ranking. Based on the above data, we can rank the performance of various DFT+vdW methods, in terms of the metrics of the simulated air−water interface. This is done by using the following equation  The ranking for g OO (r) for the method j was calculated using the average ranking for individual components of g OO (r), i.e., the set {h max , r max , h min , r min }. Using the above expression, a smaller value κ i j implies better performance by method j for quantity i. For all values except density (χ ρ ref = 1.00 g/cm 3 ) and g OO (r), the data obtained from the POLI2VS MD simulation results have been taken as the reference values. Table 3 shows the performance ranking for the DFT+vdW methods used in this study. While the average nl-vdW functional score (1.03) is close to the empirical vdW (0.87), the score of the best performing nl-vdW method, i.e., revPBE-rVV10 (0.84), is much larger than the revPBE-D3(0) functional score (0.41). This illustrates that the overall performance of DFT-nl-vdW methods is either similar or worse than empirical vdW methods. nl-vdW methods enhance the bulk water density when compared to non-vdW corrected functionals, but the interfacial description of water is either similar or worse when compared to empirical vdW methods.

CONCLUSION
We tested the efficacy of two different classes of vdW models: environment-dependent empirical corrections and nonlocal vdW functionals for a description of the air−water interface. By using a scoring scheme which ranks the DFT+vdW functionals based on bulk water density, interfacial width, interfacial free O−D group fraction, interfacial free O−D lifetime decay, and oxygen−oxygen radial distribution function, we were able to study the effect of nl-vdW methods on the air− water interface. We concluded that within the margin of error both nl-vdW methods (optB88-DRSLL and rVV10) show a similar or worse description for interfacial water as compared to the empirical D2 and D3 vdW methods. Our data clearly indicates that the revPBE+D3 models are much better suited for simulating the air−water interface, while the nl-vdW methods perform similarly to the PBE functional.
Finally, we would like to note that our results do not indicate that the nl-vdW corrections have no opportunity to describe the water dynamics accurately. Rather, our data illustrate that the current nl-vdW correction schemes are not well optimized, and thus further fine-tuning of the nl-vdW schemes for an accurate description of the air−water interface will be required.

Journal of Chemical Theory and Computation
Article