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 - 602 ] 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, Etot, within the GGA[17−19] can be expressed aswhere n() denotes the electron density, Ts[n()] denotes the kinetic energy, Eext[n()] denotes the electron-nuclei interaction energy, EH[n()] denotes the electron–electron interaction energy, EXGGA[n()] denotes the exchange energy of the GGA functional, ECGGA[n()] denotes the correlation energy of the GGA functional, and Eempirical-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 ECGGA/LDA[n()] is the correlation energy calculated within the GGA or the local density approximation (LDA), and ECnl-vdW[n()] 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 .[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 sum-frequency 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

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 aswhere C6 is the coefficient for atoms A and B, r is the distance between A and B, s6 is a scaling factor, and fdamp(r) is the short-range damping function. C6 is empirically determined and is independent of the local environment, making this correction scheme categorized into the environment-independent empirical vdW correction scheme.[52]

Environment-Dependent Empirical vdW Corrections

Straightforward extension of the Grimme’s D2 correction is to make the coefficient C dependent on the local environment. The vdW correction schemes using the environment-dependent C(CN,CN) coefficients are called the environment-dependent empirical vdW correction,[52] where CN and CN 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(CN,CN) 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]

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 ECnl-vdW[n()] can be written aswhere the correlation kernel ϕ is a function of the density n() and the gradient of the density ∇n(). 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(),|∇n()|) 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() and ∇n(). In the rVV10 approach ϕ contains two parameters - b 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 (ECLDA), while the rVV10 method uses the correlation energy calculated within the GGA (ECGGA).

Simulation Protocols

FPMD Simulation

We have performed Born–Oppenheimer simulations using the QUICKSTEP[58] method implemented in the CP2K code.[59] We have used Perdew–Burke–Ernzerhof (PBE)[60] and revised PBE (revPBE)[60,61] for the GGA XC functionals and Vosko-Wilk-Nusiar (VWN)[62] for the LDA correlation functional of the DRSLL method. This XC functional energy is combined with empirical vdW correction schemes using Grimme’s D2,[12] D3(0),[13] D3(BJ),[14] and D3m(BJ),[15] as well as the nl-vdW correction schemes using the DRSLL[20,21] and rVV10.[22,23] Furthermore, we tested the optB88-DRSLL method.[19,29] The resulting calculation methods and their notations are listed in Table .

Table 1

List of Simulations Methods Discussed

notation	exchange	correlation	vdW correction	remarks
PBEa	PBE	PBE
PBE-D2a	PBE	PBE	D2b
PBE-D3(0)a	PBE	PBE	D3(0)c
PBE-D3(BJ)a	PBE	PBE	D3(BJ)d
PBE-D3m(BJ)	PBE	PBE	D3m(BJ)e
PBE-DRSLL	PBE	VWN	DRSLLf
PBE-rVV10	PBE	PBE	rVV10g	b = 6.6, C = 0.0093
revPBEa	revPBE	PBE
revPBE-D2	revPBE	PBE	D2b
revPBE-D3(0)a	revPBE	PBE	D3(0)c
revPBE-D3(BJ)a	revPBE	PBE	D3(BJ)d
revPBE-DRSLL	revPBE	VWN	DRSLLf
revPBE-rVV10	revPBE	PBE	rVV10g	b = 3.7, C = 0.0093
optB88-DRSLL	optB88h	VWN	DRSLLf

The data are obtained from our previous work.[11]

Reference (12).

Reference (13).

Reference (14).

Reference (15).

Reference (20).

Reference (22).

Reference (29).

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 (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, D2O was used instead of H2O, 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 D2O molecules in a simulation cell (L,L,L) = (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.

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−74] 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

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)where zG 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 D2O 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/cm3)[11] is very close to the density predicted in the NPT ensemble (0.993 ± 0.001 g/cm3)[66] in the MD simulation with the POLI2VS force field model.

Fraction of Interfacial Water Molecules with the Free O–D Group and the Angle of the Free O–D Group

An O–D group of the D2O 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 z-axis 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,75] 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,47−49,76−78]

Lifetime of the Free O–D Group

The lifetime of the free O–D group is computed using the time correlation functionwhere h(t) is defined as 1 when an O–D group is free at time t, otherwise 0. C(t) is fitted using the double exponential form[79]within the time region 0–5 ps, where a, b, and c are the fitting coefficients. The time constants τf and τs denote the fast and slow decays, respectively. τf is governed by the librational motion of water, while τs is governed by the motion that the free O–D group rotates and forms a hydrogen bond at the water–air interface. These values can be probed using time-resolved SFG spectroscopy,[41,73,80] whereas so far only the lifetime of the free O–H group has been obtained[41] and the lifetime of the free O–D group is not determined. As the POLI2VS simulation agrees with the experimental data for H2O,[32] we used the POLI2VS data of D2O as the reference.

Radial Distribution Function (RDF)

The oxygen–oxygen RDF gOO(r) was calculated when the reference oxygen atom of a water molecule is in the region |z| ≤ 1.35 Å of the slab model. The RDF was calculated based on the computed ρ0. The resolution of the RDF, Δr, was set as 0.1 Å.

Results

The obtained data is summarized in Table with the earlier evaluated values for the PBE and revPBE with and without vdW corrections.

Table 2

Bulk and Interfacial Water Data Using Various DFT+vdW Methodsa

		PBE							revPBE
	optB88-DRSLL	none	D2	D3(0)	D3(BJ)	D3m(BJ)	DRSLL	rVV10	none	D2	D3(0)	D3(BJ)	DRSLL	rVV10	POLI2VS	exp.
ρ (g/cm3)	1.07	0.91	1.02	1.03	1.00	1.00	1.12	1.03	0.69	0.94	0.94	0.93	0.93	1.06	0.99	1.00
zG (Å)	8.08	9.50	8.43	8.44	8.70	8.63	7.70	8.36	12.5	9.22	9.18	9.31	9.31	8.14	8.76
δ (Å)	1.02	1.08	1.06	1.01	1.05	1.02	0.81	0.96	2.06	1.06	1.17	1.16	1.06	1.03	1.29
DA (%)	7	12	7	8	9	9	6	7	27	12	11	14	19	9	12
DAA (%)	13	15	12	13	14	13	14	13	11	14	13	14	17	13	16
DDA (%)	16	17	17	17	17	17	17	17	13	18	17	18	19	17	17
free OH (%)	20	27	19	21	22	22	20	19	38	26	24	28	35	23	28
<ϑ> (deg)	52	54	51	53	51	55	56	48	76	65	66	65	72	56	61
τf (ps)	0.15	0.19	0.16	0.15	0.18	0.17	0.09	0.21	0.19	0.12	0.11	0.15	0.10	0.14	0.11
τs (ps)	1.69	2.18	1.99	1.87	2.04	1.86	1.00	2.42	1.05	1.03	1.06	1.25	0.51	1.48	1.02	1.1b
hmax	3.12	3.67	3.62	3.44	3.78	3.69	2.54	3.73	2.74	2.97	2.72	3.12	2.18	2.95	2.73	2.58c
rmax (Å)	2.75	2.75	2.71	2.72	2.71	2.71	2.82	2.70	2.86	2.77	2.81	2.77	2.92	2.74	2.78	2.80c
hmin	0.62	0.18	0.31	0.36	0.18	0.21	−	0.23	0.46	0.56	0.69	0.42	−	0.67	0.83	0.84c
rmin (Å)	3.24	3.36	3.25	3.26	3.27	3.28	−	3.24	3.61	3.42	3.46	3.36	−	3.23	3.43	3.47c

The average error bars for ρ, zG, δ, fractions, ⟨ϑ⟩, τf, τs, hmax, rmax, hmin, and rmin are 0.01 g/cm3, 0.1 Å, 0.2 Å, 1.0%, 1.5 degrees, 0.01 ps, 0.2 ps, 0.1, 0.01 Å, 0.4, and 0.02 Å, respectively.

Reference (32).

Reference (84).

Bulk Water Density

First, we compare the obtained density values with the reported value. As is mentioned in ref (11), 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/cm3) is similar to the reported values of 1.13 g/cm3 in ref (26) and 1.18 g/cm3 in ref (28). However, our data for the revPBE-DRSLL differs from the data in ref (26); we obtained ρ = 0.93 ± 0.01 g/cm3, smaller than the value of 1.02 g/cm3.[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/cm3. 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 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/cm3.

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.

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 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 . Figure displays the time correlation function C(t) for the DFT XC functionals, and Table summarizes the fitted time constants. The fast time scale (τf) is in the range of 0.1–0.2 ps, while the slow time scale (τs), which is associated with the motion of the free O–D rotation and formation of a new hydrogen bond, is scattered from 0.5 to 2.4 ps. We discuss the behavior of the slow time scale in more detail.

Figure 1

Time correlation functions of the free O–D group of the interfacial water molecules.

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.[81,82] 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.

RDF

We plot gOO(r) in Figure and summarize the heights and positions of the first maxima (hmax;rmax) and minima (hmin;rmin) in Table . 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,83] 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 hmax = 2.77 and rmax = 2.80 Å at 295 K in ref (24) and hmax = 2.51 and rmax = 2.80 Å at 295 K in ref (28). These agree with hmax = 2.54 and rmax = 2.82 Å in our simulation. hmax = 2.18 and rmax = 2.92 Å in our simulation at the revPBE-DRSLL level of theory are corroborated by previously reported values of hmax = 2.35 and rmax = 2.92 Å at 300 K in ref (24). Furthermore, one can see that gOO(r) simulated with the DRSLL correction scheme does not show a clear RDF minimum, indicating that the water structure is disordered excessively. This is consistent with ref (26).

Figure 2

Oxygen–oxygen RDF of water simulated in the FPMD simulation with the (a) PBE XC functional, (b) revPBE XC functional, and (c) optB88-DRSLL level of theory. The shaded area represents the experimental data from X-ray diffraction measurements.[84]

The data in Table 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 hmax.

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/cm3 at 295 K and 1 atm,[25] which is consistent with our density value of ρ = 1.07 ± 0.02 g/cm3. This density value is higher than the experimental value of ρ = 1.00 g/cm3. 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 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.

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 equationwhere κ is the score of the method j = {optB88, revPBE, revPBE-D3(0), revPBE-D3(BJ), revPBE-DRSLL, revPBE-rVV10, PBE, PBE-D3(0), PBE-D3(BJ), PBE-D3m(BJ), PBE-DRSLL, PBE-rVV10} for the calculated property i = {ρ, interfacial free O–D fraction, ⟨θ⟩, τs, gOO(r)}. χ 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 of quantity i for all methods in set j. The ranking for gOO(r) for the method j was calculated using the average ranking for individual components of gOO(r), i.e., the set {hmax, rmax, hmin, rmin}. Using the above expression, a smaller value κ implies better performance by method j for quantity i. For all values except density (χρref = 1.00 g/cm3) and gOO(r), the data obtained from the POLI2VS MD simulation results have been taken as the reference values.

Table 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.

Table 3

Score of Various DFT+vdW Functionals for Bulk and Interfacial Data

		PBE							revPBE
	optB88-DRSLL	none	D2	D3(0)	D3(BJ)	D3m(BJ)	DRSLL	rVV10	none	D2	D3(0)	D3(BJ)	DRSLL	rVV10
ρ	0.66	0.90	0.22	0.28	0.05	0.02	1.16	0.32	2.98	0.59	0.57	0.70	0.70	0.60
δ	0.96	0.73	0.81	0.96	0.83	0.96	1.67	1.16	2.69	0.80	0.40	0.47	0.82	0.91
fraction	1.33	0.15	1.46	1.22	0.94	1.04	1.40	1.44	1.75	0.29	0.62	0.04	1.26	0.90
<ϑ>	1.05	0.85	1.27	0.95	1.26	0.81	0.58	1.52	1.71	0.39	0.47	0.37	1.25	0.62
τs	1.20	2.09	1.74	1.52	1.83	1.51	0.04	2.51	0.06	0.02	0.07	0.42	0.91	0.83
hmax	1.08	2.15	2.07	1.71	2.38	2.21	0.07	2.27	0.32	0.79	0.28	1.08	0.78	0.75
rmax	0.77	0.77	1.39	1.24	1.39	1.39	0.31	1.54	0.93	0.46	0.15	0.46	1.85	0.93
hmin	1.16	3.43	2.75	2.49	3.41	3.28	−	3.17	1.97	1.45	0.75	2.17	−	0.88
rmin	1.98	0.95	1.89	1.81	1.72	1.63	−	1.98	1.20	0.43	0.09	0.95	−	2.06
final score	1.07	1.09	1.25	1.12	1.19	1.08	0.82	1.53	1.72	0.48	0.41	0.53	0.93	0.84

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.