Two hybrid van der Waals density functionals (vdW-DFs) are developed using 25% Fock exchange with (i) the consistent-exchange vdW-DF-cx functional [K. Berland and P. Hyldgaard, Phys. Rev. B 89, 035412 (2014)] and (ii) with the vdW-DF2 functional [K. Lee et al., Phys. Rev. B 82, 081101 (2010)]. The ability to describe covalent and non-covalent binding properties of molecules is assessed. For properties related to covalent binding, atomization energies (G2-1 set), molecular reaction energies (G2RC set), and ionization energies (G21IP set) are benchmarked against experimental reference values. We find that hybrid-vdW-DF-cx yields results that are rather similar to those of the standard non-empirical hybrid PBE0 [C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999)], with mean average deviations (MADs) of 4.9 and 5.0 kcal/mol for the G2-1 set, respectively. In this comparison, experimental reference values are used, back corrected by wavefunction-based quantum-chemistry calculations of zero-point energies. Hybrid vdW-DF2 follows somewhat different trends, showing on average significantly larger deviations from the reference energies, with a MAD of 14.5 kcal/mol for the G2-1 set. Non-covalent binding properties of molecules are assessed using the S22 benchmark set of non-covalently bonded dimers and the X40 set of dimers of small halogenated molecules, using wavefunction-based quantum chemistry results as references. For the S22 set, hybrid-vdW-DF-cx performs better than standard vdW-DF-cx for the mostly hydrogen-bonded systems, with MAD dropping from 0.6 to 0.3 kcal/mol, but worse for purely dispersion-bonded systems, with MAD increasing from 0.2 to 0.6 kcal/mol. Hybrid-vdW-DF2 offers a slight improvement over standard vdW-DF2. Similar trends are found for the X40 set, with hybrid-vdW-DF-cx performing particularly well for binding energies involving the strongly polar hydrogen halides, but poorly for systems with tiny binding energies. Our study of the X40 set reveals the potential of mixing Fock exchange with vdW-DF, but also highlights shortcomings of the hybrids constructed here. The solid performance of hybrid-vdW-DF-cx for covalent-bonded systems, as well as the strengths and issues uncovered for non-covalently bonded systems, makes this study a good starting point for developing even more accurate hybrid vdW-DFs.

Hybrid functionals are widely used for computing molecular properties because they generally offer more accurate thermochemical and structural properties than non-hybrid semilocal functionals in density functional theory (DFT).1–8 Through various procedures, hybrid functionals mix in a fraction of Fock exchange, introducing a Fock term in the exchange part of a Kohn-Sham density functional. Many semi-empirical or fitted hybrid functionals have been developed, with B3LYP being one of the most popular.1,2,8,9 One of the most successful hybrids, the PBE0 functional,3 on the other hand, mixes in 25% of Fock exchange with the exchange of the PBE,10,11 a fraction which is motivated by ongoing arguments from many-particle theory.12–14 

In many molecular complexes, London dispersion forces contribute significantly to the binding energy of the system. To accurately describe all binding properties of such a system, it is desirable to include both a non-local correlation and a Fock-exchange fraction within the same functional.8 Non-local correlations are missing in both semi-local functionals in the generalized-gradient approximation (GGA) and in standard hybrid functionals, but they are included from the onset in the van der Waals density functional (vdW-DF) method.15–22 vdW-DFs have found a wide user base in physics and material science.21,23,24 They are regularly used in the study of adsorption,25–29 layered systems,30–34 molecular crystals,35–37 and more recently solids with covalent and metallic bonds.20,38,39 However, assessments of how well vdW-DFs describe general properties of molecular systems are scarce, which is likely related to the fact that hybrid vdW-DFs have so far gone unexplored.

In this paper, we define and assess the performance of two simple, unscreened hybrid van der Waals functionals, termed vdW-DF-cx0 and vdW-DF2-0, for molecular systems. A number of dispersion-corrected (DFT-D) hybrid functionals do already exist.40–43 Those descriptions rely on atom-centered pair-potential corrections to account the non-local correlation. In more advanced variants, the coefficients of these pair potentials can depend on, for instance, the electron density, as they are in DFT-TS and more recent extensions44,45 and in associated hybrids.46 Hybrid variants of the Vydrov-Voorhis (VV) nonlocal density functional,47 which borrows several features of Tkatchenko-Scheffler formulation of vdW-corrected DFT, vdW-TS, has also been developed. All of these schemes include one or several adjustable parameters that can be re-tuned when combined with a hybrid functional to achieve good performance for a set of reference systems. Here, we investigate the performance of two hybrid vdW-DFs where we deliberately avoid fitting to reference systems.

The vdW-DF framework is formally well anchored in many-body physics.15,16,19,21,48 Different versions of the functional have been developed,15,38,49,50 such as vdW-DF217 and vdW-DF-cx,18,24 which emphasize different limits of many-body theory. In this work, we define two new vdW-DF hybrids, vdW-DF-cx0 and vdW-DF2-0, modifying vdW-DF-cx and vdW-DF2 to introduce 25% of Fock exchange, in analogy with the design of PBE0.3,12–14 The 25% mixing fraction has been rationalized1,12–14 in terms of a coupling-constant λ integration in the adiabatic-connection formula (ACF).51,52 In a mean-value evaluation of the ACF, using a fraction of Fock exchange takes care of the λ0 limit, where the functional is exclusively represented by exchange. At the λ1 end, the properties of a system are dominated by plasmons,53 suggesting that functionals of the electron-gas tradition, like PBE, provide robust starting points for hybrid constructions, like PBE0.12,13 However, the same rationale also applies to the consistent-exchange vdW-DF-cx version, which rely on vdW-DF’s single plasmon-pole description for specifying the exchange and correlation, except in very inhomogeneous regions.18,24 In fact, vdW-DF, and this version in particular, can be viewed as a series expansion of the effective response defined by a GGA-type exchange-correlation hole.19,21

Even if less strictly enforced, vdW-DF2 is also designed with consistency between the exchange and correlation in mind and can thus also be seen as plasmon based.17,19,21,24 Whereas vdW-DF-cx has a good consistency between its rather soft exchange functional and the internal exchange that sets the strength of the non-local correlation term, vdW-DF2 relies on more rapidly increasing exchange enhancement factors both in the exchange and within the non-local correlation.54,55 vdW-DF-cx and vdW-DF2 can therefore be viewed as two different limits of viable vdW-inclusive density functionals. Thus by benchmarking these two functionals for molecular systems, we gain insight into how to best construct an accurate hybrid vdW-DF.

vdW-DF-cx0 and vdW-DF2-0 are implemented within version 6.0 of the quantumespresso package.56 This version includes an implementation of spin-vdW-DF,20 which is needed to correctly compute atomization energies. In all calculations, Martins-Troullier57 norm-conserving PBE pseudopotentials of the FHI-abinit58 kind are employed. We use an energy cutoff of 80 Ry. Atomic coordinates were relaxed (to 0.0025 eV/Å) in all studies, except for the S22 set.59 Like in several earlier studies,17,18 the internal molecular coordinates are kept fixed for the S22 set and the optimal separations are determined by calculating potential energy curves.

A good hybrid vdW-DF should be able to describe both covalent and non-covalent binding accurately. To assess the ability of vdW-DF-cx0 and vdW-DF2-0 to describe covalent binding properties, we consider atomization energies (in the G2-1 subset) of molecules and chemical reaction energies (in the G2RC subset) of the G2 benchmark set.60–62 We furthermore benchmark molecular ionization potentials (the G2IP subset).60 We compare with experimental reference values, in the case of G2-1 and G2RC subsets back-corrected by wavefunction-based quantum chemistry results for zero-point energies.60 

Figure 1 provides an overview of the mean absolute deviation (MAD) of the energies of the three sets, displaying results of vdW-DF-cx0, vdW-DF2-0, and their non-hybrid equivalents.62 For reference, PBE and PBE0 results are also shown. The comparison shows that among the non-hybrid functionals, PBE performs best, except for atomization energies where vdW-DF2 performs better. Among the hybrids, vdW-DF-cx0 has the best overall performance, slightly better than that of PBE0. Both vdW-DF-cx0 and PBE0 perform better than the corresponding non-hybrid functionals. In contrast, vdW-DF2-0 does not necessarily perform better than vdW-DF2. To gain insight into these trends, we next consider the results of the three sets in more detail.

FIG. 1.

Mean absolute deviations of the three benchmark sets considered for covalent molecular properties.

FIG. 1.

Mean absolute deviations of the three benchmark sets considered for covalent molecular properties.

Close modal

Table I gives an overview of our results for atomization energies, in terms of MAD, mean absolute relative deviation (MARD), mean deviation (MD), and mean relative deviation (MRD). It shows that vdW-DF2-0 is significantly less accurate than vdW-DF2, but vdW-DF-cx0 is significantly more accurate than vdW-DF-cx. Comparing the MRD values of a non-hybrid with a corresponding hybrid indicates that the atomization energies are on average about 10% lower for the hybrid than for the non-hybrids. Since vdW-DF2 underestimates atomization energies on average by about 2%, vdW-DF2-0 underestimates them by about 12%. This systematic underestimation leads to large MAD and MARD values of vdW-DF2-0. vdW-DF-cx overestimates atomization energies by 9% on average, giving it the worst performance among the non-hybrid functionals. However, this overestimation is a good match with the reduction in atomization energies when introducing 25% Fock exchange, resulting in an MRD value of merely −1.5%. In turn, vdW-DF-cx0 ends up with the smallest MARD value (3.7%) among the tested functionals, somewhat smaller than the PBE0 value 4.1%.

TABLE I.

Benchmark of the atomization energies in the G2-1 set.60 The table lists mean relative deviation (MRD), mean absolute relative deviation (MARD), mean deviation (MD), and MAD values of the functionals considered here relative to the experimental references.

MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) 3.6 −2.6 −0.8 −7.2 5.8 −1.4 
MARD (%) 5.9 4.1 4.6 9.5 7.4 3.7 
MD (kcal/mol) 5.8 −3.4 −2.3 −12.1 9.1 −1.5 
MAD (kcal/mol) 7.5 5.0 6.6 14.5 9.7 4.9 
MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) 3.6 −2.6 −0.8 −7.2 5.8 −1.4 
MARD (%) 5.9 4.1 4.6 9.5 7.4 3.7 
MD (kcal/mol) 5.8 −3.4 −2.3 −12.1 9.1 −1.5 
MAD (kcal/mol) 7.5 5.0 6.6 14.5 9.7 4.9 

For reaction energies, the hybrid vdW-DFs outperform the non-hybrids. While Table II shows that the MAD becomes somewhat smaller, the relative deviations are reduced significantly. For instance, the MARD value of 47% for vdW-DF-cx drops to 30% for vdW-DF-cx0. For vdW-DF2 and vdW-DF2-0, the corresponding numbers are 59% and 36%. Despite having poorer atomization energies, as discussed earlier, vdW-DF2-0 predicts reaction energies more accurately than vdW-DF2. This suggests a partial error cancellation between the atomization energies of both the product and the reactant.

TABLE II.

Comparison of the performance of PBE0 and the hybrid van der Waals density functionals, vdW-DF2-0 and vdW-DF-cx0, for chemical reaction energies in the G2RC set.60 

MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) −15.8 5.3 −31.1 −11.4 −26.4 −8.6 
MARD (%) 37.3 33.0 59.4 36.0 47.3 30.2 
MD (kcal/mol) 0.7 −2.8 6.8 2.7 1.9 −1.1 
MAD (kcal/mol) 5.6 6.1 10.0 5.8 6.1 4.8 
MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) −15.8 5.3 −31.1 −11.4 −26.4 −8.6 
MARD (%) 37.3 33.0 59.4 36.0 47.3 30.2 
MD (kcal/mol) 0.7 −2.8 6.8 2.7 1.9 −1.1 
MAD (kcal/mol) 5.6 6.1 10.0 5.8 6.1 4.8 

Figure 2 shows that the deviations grow quite slowly with the size of the reaction energies. In the figure, the reactions are ordered according to the size of their respective reaction energies (gray bars in the upper and middle panels, magnitude given by the right vertical axis). The comparison shows that the deviations of vdW-DF-cx0 are typically quite similar to those of PBE0, while overall being slightly more accurate. vdW-DF-cx also exhibits deviations that are fairly similar to those of PBE (as shown in the top panel of Fig. 2). In contrast, the deviations of vdW-DF2 and vdW-DF2-0 do not follow the same general trends as PBE, PBE0, vdW-DF-cx, and vdW-DF-cx0.

FIG. 2.

Reaction energies of the G2RC subset of the G2 sets60 for the different functionals compared here, as ordered by reference reaction energies. These energies are defined as the energy difference between reactants and products and are indicated by the gray bars in the background (right vertical axis), on a scale five times larger than the deviations. We report deviations, that is, the differences between our calculated results and the reference data. The dashed gray line is a guide to the eye for the magnitude of the reaction energy. The lowest panel shows the corresponding relative deviations, where the inset provides a smaller span on the vertical axis, for the systems with reference values larger than 20 kcal/mol.

FIG. 2.

Reaction energies of the G2RC subset of the G2 sets60 for the different functionals compared here, as ordered by reference reaction energies. These energies are defined as the energy difference between reactants and products and are indicated by the gray bars in the background (right vertical axis), on a scale five times larger than the deviations. We report deviations, that is, the differences between our calculated results and the reference data. The dashed gray line is a guide to the eye for the magnitude of the reaction energy. The lowest panel shows the corresponding relative deviations, where the inset provides a smaller span on the vertical axis, for the systems with reference values larger than 20 kcal/mol.

Close modal

The bottom panel of Fig. 2 shows the relative deviations of the different reaction energies. The relative deviations are typically less than 20% for the reactions with large reaction energies, as shown in the inset of the bottom panel, but they can be much larger for the cases with small reaction energies. In the worst cases, even the sign can be wrong, corresponding to relative deviations below −100%. The hybrid variants generally reduce the deviations of the reaction energies. As an example, the energy of the reaction CH2O2CO2+H2 (shown in the inset) is −2.0 kcal/mol. vdW-DF2 incorrectly predicts this to be an endothermic reaction with a reaction energy of 6.1 kcal/mol which drops to 2.9 kcal/mol when going to vdW-DF2-0. vdW-DF-cx correctly predicts this reaction to be exothermic, but only barely so, with a reaction energy of −0.1 kcal/mol, and vdW-DF-cx0 however gives a value of −2.3 kcal/mol which is only a 12% overestimation of the reference energy. The improvements for the systems with small reaction energies are the main reason why the MARD value is significantly lower for the hybrids, while the reduction in MAD is more modest.

For ionization potentials, the comparison for the G21IP set, summarized in Table III and Fig. 1, shows that both the hybrids perform similar or slightly better. vdW-DF-cx and vdW-DF-cx0 perform better than vdW-DF2 and vdW-DF2-0.

TABLE III.

Comparison of the performance of PBE0 and the hybrid van der Waals density functionals, vdW-DF2-0 and vdW-DF-cx0, for ionization potential values in the G21IP set.60 

MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) −0.4 −0.2 3.1 3.0 −0.6 0.2 
MARD (%) 2.0 2.0 3.3 3.3 1.7 1.7 
MD (kcal/mol) −1.3 −0.9 7.4 7.1 −1.8 0.1 
MAD (kcal/mol) 4.9 4.7 7.8 7.7 4.2 3.9 
MoleculePBEPBE0DF2DF2-0DF-cxDF-cx0
MRD (%) −0.4 −0.2 3.1 3.0 −0.6 0.2 
MARD (%) 2.0 2.0 3.3 3.3 1.7 1.7 
MD (kcal/mol) −1.3 −0.9 7.4 7.1 −1.8 0.1 
MAD (kcal/mol) 4.9 4.7 7.8 7.7 4.2 3.9 

Compared to other methods, vdW-DF-cx0 performs quite well for G2RC and G21IP. The three hybrids, B3LYP, B3PW91, and revPBE0 amended with VV10 correlation,63 give MADs of, respectively, 4.3, 8.3, and 8.5 kcal/mol for G2RC and 5.2, 4.5, and 6.0 kcal/mol for G2IP.

The performance of vdW-DF-cx0 for covalent molecular binding properties makes this a promising vdW-DF hybrid functional candidate. vdW-DF2-0 is not as promising, which is clearly illustrated by its atomization energies. In a separate test, we find that the difference between the trends of vdW-DF2-0 and vdW-DF-cx0 for small-molecule atomization energies is predominantly due to the exchange component. For covalent interactions, the performance is rather insensitive to whether vdW-DF1 or vdW-DF2 correlation is employed for a fixed exchange functional. This is in contrast to the case for non-covalent interactions.54 

Figure 3 shows the results for the S22 set59 following the standard division into hydrogen-bonded, dispersion-bonded, and mixed systems. vdW-DF2-0 has a smaller MAD than vdW-DF2 for all the three cases. vdW-DF-cx0 and vdW-DF-cx have a similar MAD of about 0.4 kcal/mol. Most striking, however, is the poorer performance of vdW-DF-cx0 compared to vdW-DF-cx for purely dispersion-bonded systems—with MAD increasing from 0.2 to 0.6 kcal/mol—combined with an improved performance for the hydrogen bonded systems, MAD decreasing from 0.6 to 0.3 kcal/mol. The poor results for dispersion bonded systems are due to overestimation of binding energies.

FIG. 3.

The performance of the functionals considered here for the S22 set non-covalent.

FIG. 3.

The performance of the functionals considered here for the S22 set non-covalent.

Close modal

The improved performance for hydrogen bonded systems is linked to the fraction of Fock exchange, as the performance also improves correspondingly for vdW-DF2 and PBE. The MAD shrinks from 1.2 kcal/mol to 1.1 kcal/mol when going from PBE to PBE0. However, these numbers are inferior to the vdW-DFs as dispersion forces generally contribute to the binding in all kinds of non-covalently bonded dimers, including systems that are mostly hydrogen bonded. We therefore omit PBE and PBE0 in further analysis of non-covalently bonded systems.

In the community there are two typical ways of assessing the accuracy of the S22 dataset: one is to rely on coordinates fixed to the optimum coordinates of the benchmark set59 and the other is to relax inter-atomic positions. While both have pros and cons, we prefer the latter, as using static coordinates introduces a slight preference for functionals with overestimated binding energies. As long as the optimum separation does not exactly agree with the reference geometry, one would always compare with an energy that is smaller than the optimum one. For reference with literature, the MAD values calculated at fixed reference geometries are 0.7 and 0.4 kcal/mol for vdW-DF2 and vdW-DF2-0 and 0.4 and 0.3 kcal/mol for vdW-cx and vdW-DF-cx0, respectively.

In assessing the accuracy of the different methods, we need to keep in mind that our calculations are based on the non-native norm-conserving PBE pseudopotential. This is a common practice in vdW-DF calculation when using plane-wave codes such as Quantum Espresso, with good cross-code transferability,64 so it gives a good indication of the accuracy obtained in practice. However, native pseudopotentials or all-electron calculations would be preferable. To get an indication of the inaccuracy introduced by using norm-conserving PBE pseudopotentials, our vdW-DF2 results with geometries of the S22 set can be compared with corresponding results generated by Vydrov and van Voorhis47 using a Gaussian basis at the aug-cc-pVTZ level. Relative to this reference, our vdW-DF2 results have a mean absolute error of 0.2 kcal/mol and mean error of −0.1 kcal/mol (mean absolute relative error of 2%). As the mean absolute error is comparable to the accuracy of the hybrid vdW-DF results for the S22 set, the precision is limited. In gauging the accuracy of the MAD, however, the magnitude of the mean error of 0.1 kcal/mol gives a better indication of the inaccuracy of the MAD as the two sources of deviations are not related.

Vydrov and van Voorhis47 also showed that the long-range corrected (LC) Fock exchange variant of VV10 also exhibits improved performance for hydrogen bonding compared to VV10, with MAD dropping from 0.7 kcal/mol to 0.3 kcal/mol, while for mixed and dispersion bonded system, standard VV10 is in fact slightly more accurate. With an overall MAD of 0.2 kcal/mol for the S22 set, LC-VV10 is more accurate than the vdW-DF hybrids considered here, which is not surprising given that VV10 and LC-VV10 are reference systems optimized to S22 binding curves. Other variants of VV10, combined with B3LYP, B3PW91, and revPBE0 exchange-correlation, result in MAD values of 0.5, 0.5, and 0.2 kcal/mol, respectively.63 

The clear overestimation of vdW-DF-cx0 for dispersion bonded systems shows that Fock exchange and vdW-DF-cx exchange behave quite differently for dispersion-bonded systems, as vdW-DF-cx itself has a rather low MAD. Thus the vdW-DF1 correlation (which is also used in vdW-DF-cx) may not be the optimal one for the exchange description used in vdW-DF-cx0. The low MAD for hydrogen-bonded systems nonetheless indicates the promise of mixing a soft exchange like that of vdW-DF-cx with Fock exchange.

The mixed vdW-DF-cx0 performance for the S22 set motivates us to also consider the X40 set65 of non-covalently bonded dimers involving at least one halogenated molecule. This set of dimers involves small molecules like methane, and its halogenated derivatives, an example being the methane-Cl2 dimer (reference binding energy of 1.08 kcal/mol) or bromo-benzene-acetone dimer (2.43 kcal/mol). They also include the strongly polar hydrogen halides, such as the HF-methanol dimer (9.59 kcal/mol).

Figure 4 shows the deviation in binding energies (upper panel) and binding separation (lower panel) ordered by the magnitude of the reference energy. The bars in the upper panel indicate the magnitude of the reference energies, with scale on the right vertical axis. The bar color indicates a sorting of the dimers into molecules involving methane (pink), benzene (green), hydrogen halide (orange), or none of these (gray). We organize the X40 set into these subsystems for the following reasons: For all the dimers of methane and a “halogenated molecule,” the binding energies are significantly overestimated so it overshadows all other MRD trends. For benzene-“halogenated molecule” dimers, the aromatic π-orbitals could contribute to the binding, giving them a somewhat different binding nature. For the dimers involving hydrogen halides, the binding is dominated by the strong polarity of the halide.

FIG. 4.

Our benchmark of the X40 set of molecular dimers involving halogenated molecules ordered by the reference energies. The connected markers in the upper panel indicate the deviations in the binding energy of the various molecular dimers, that is, the differences between our calculated results and the reference values. The bars show the reference energy on a ten times larger scale (right vertical axis). The dimers have been sorted into systems involving methane (pink bars), involving benzene (green bars), involving hydrogen halides (orange), and the remainders (gray). The dashed curves correspond to +/−, respectively, 10% and 30% of the binding energy. The lower panel shows the performance for intermolecular separations, defined in terms of the shortest atomic separations between dimers. The bars indicate the full reference separations on a ten times larger scale (right vertical axis).

FIG. 4.

Our benchmark of the X40 set of molecular dimers involving halogenated molecules ordered by the reference energies. The connected markers in the upper panel indicate the deviations in the binding energy of the various molecular dimers, that is, the differences between our calculated results and the reference values. The bars show the reference energy on a ten times larger scale (right vertical axis). The dimers have been sorted into systems involving methane (pink bars), involving benzene (green bars), involving hydrogen halides (orange), and the remainders (gray). The dashed curves correspond to +/−, respectively, 10% and 30% of the binding energy. The lower panel shows the performance for intermolecular separations, defined in terms of the shortest atomic separations between dimers. The bars indicate the full reference separations on a ten times larger scale (right vertical axis).

Close modal

We also show values for unrelaxed coordinates for reference with literature.

Table IV summarizes the trends for the X40 set and for various subsystems. For the dimers involving methane, all the functionals overestimate binding energies, with overall somewhat better performance for the hybrids. The dimers involving methane in the X40 set are always paired with another small molecule, either a diatomic halogen going from F2 to I2 (system 1, 5, 8, and 10) or halogenated methane derivatives (system 2, 3, 4, and 6). The overestimation is evident both in the table and by cases in which the data are outside either the first or second dashed curve indicating 10% and 30% deviation in the upper panel of Fig. 4. vdW-DF-cx and vdW-DF-cx0 also significantly overestimate binding separations of the methane systems, on average by, respectively, 0.25 Å and 0.17 Å, as with vdW-DF-cx for the methane dimer of the S22 set.18 Issues for the smallest molecules like methane is a known shortcoming of vdW-DF-cx, which is a consequence of using the Langreth-Vosko exchange66 form in a wide regime of reduced gradients s, making it overly repulsive for systems dominated by large s values.18 The separations of vdW-DF2 and vdW-DF2-0 are rather similar for all the systems and both are on average overestimated by about 0.05 Å. This result can be traced to the fact that for molecular dimers, PW86r exchange and Hartree-Fock exchange give quite similar binding curves.67 

TABLE IV.

Summary of our benchmark of the X40 set binding energies.

DF2DF2-0DF-cxDF-cx0
Relaxed coordinates     
Involving methane     
MAD (kcal/mol) 0.3 0.2 0.2 0.2 
MARD (%) 39 20 26 30 
Involving benzene     
MAD (kcal/mol) 0.1 0.1 0.3 0.5 
MARD (%) 10 
Involving hydrogen halides     
MAD (kcal/mol) 0.3 0.6 0.5 0.1 
MARD (%) 
Other dimers     
MAD (kcal/mol) 0.2 0.2 0.3 0.2 
MARD (%) 10 10 
Full set     
MAD (kcal/mol) 0.2 0.2 0.3 0.2 
MARD (%) 13 12 12 
Unrelaxed coordinates     
MAD (kcal/mol) 0.2 0.2 0.3 0.3 
MARD (%) 12 10 
DF2DF2-0DF-cxDF-cx0
Relaxed coordinates     
Involving methane     
MAD (kcal/mol) 0.3 0.2 0.2 0.2 
MARD (%) 39 20 26 30 
Involving benzene     
MAD (kcal/mol) 0.1 0.1 0.3 0.5 
MARD (%) 10 
Involving hydrogen halides     
MAD (kcal/mol) 0.3 0.6 0.5 0.1 
MARD (%) 
Other dimers     
MAD (kcal/mol) 0.2 0.2 0.3 0.2 
MARD (%) 10 10 
Full set     
MAD (kcal/mol) 0.2 0.2 0.3 0.2 
MARD (%) 13 12 12 
Unrelaxed coordinates     
MAD (kcal/mol) 0.2 0.2 0.3 0.3 
MARD (%) 12 10 

For the dimers involving benzene, vdW-DF2 and vdW-DF2-0 show excellent performance. vdW-DF-cx0 performs worst, which stems from a poor description of the trifluorobenzene-benzene dimer and the hexafluorobenzene-benzene dimer. For the dimers involving hydrogen halides, which bind strongly, vdW-DF-cx0 exhibits a MARD of only 2%. Even if the binding separations are a bit underestimated, vdW-DF-cx0 is the only functional that does not result in a large deviation of the binding energy of the HCl-methylamine dimer. This is a case in which the binding energy shrinks once Fock exchange is mixed in.

Considering the remaining 20 dimers, both hybrids exhibit a clear improvement over the non-hybrid functionals with MARD dropping by about 40% for both vdW-DF-cx0 and vdW-DF2-0. These dimers involve a range of small organic and halogenated dimers with varying binding energy and polarity. For the larger molecules, the overestimation of separations by both vdW-DF-cx and vdW-DF-cx0 turns into an underestimation in energies. vdW-DF-cx0 typically predicts a similar but somewhat more accurate separations than vdW-DF-cx.

Evaluating the entire X40 set, we find that hybrid variants do perform better than the non-hybrid functionals. Comparing the MAD and MARD values and Fig. 4 shows that vdW-DF2-0 mostly improves the performance over vdW-DF2 for the systems with small binding energy. vdW-DF-cx0 can worsen the binding energies for such systems, which is also linked to the poorer account of dispersion bonded systems for the S22 set. It can also improve the performance for systems with strong polar binding, which is linked to the good performance for hydrogen-bonded systems of the S22 set.

We have tested the performance of two hybrid vdW-DFs, assessing both covalent and non-covalent molecular binding properties. For covalent molecular binding properties, we find that vdW-DF-cx0 gives results very similar to PBE0, with slightly better performance, whereas vdW-DF2-0 is overall less accurate. This result indicates that the exchange of vdW-DF-cx0 is here a good match to the correlation component. In turn, this indicates that vdW-DF-cx can indeed, for these properties, serve the same role as PBE in mean-value evaluation of the adiabatic connection formula. Using hybrid functionals also improves the non-covalent interactions for systems where hydrogen bonding and other polar effects contribute significantly to the binding. However, vdW-DF-cx0 overestimates binding energies of purely dispersion-bonded systems of the S22 set.

Insight into a possible mechanism behind the finite vdW-DF-cx0 deviations in describing binding energies of dispersion bonded systems can be obtained by recalling the construction of vdW-DF-cx.18 The vdW-DF-cx version rests on a consistent-exchange argument whereby exchange and correlation are brought in balance by effectively using the same plasmon-pole propagator to define both terms.19,24 However, in going from vdW-DF-cx to vdW-DF-cx0, we only replaced the exchange contributing to the total energy but did not update the internal exchange within the vdW-DF-cx correlation kernel, and to some extent, we thereby broke the internal consistency of the hybrid functional. Purely dispersion-bonded systems are precisely the regime where the parameterization of the internal functional is crucial for the total binding energies of the systems.54 Further progress in constructing more accurate hybrid vdW-DF versions may therefore be achieved by ameliorating this broken consistency. Such an investigation is beyond the scope of the present paper.

Our study shows that hybrid vdW-DFs are promising general-purpose functionals for describing binding in molecules. vdW-DF-cx0 stands out as an excellent functional for describing covalent and hydrogen/halogen-bonding of molecules, but the simple mixing of 25% Fock exchange does not preserve the good performance of vdW-DF-cx for typical dispersion bonded dimers nor does it resolve issues of vdW-DF-cx for binding between the smallest molecules. This study should therefore motivate further development of hybrid vdW functionals as well as tests of other classes of materials, both dense traditional materials and sparse vdW materials. In particular, materials, where interactions compete and both the sparse- and dense-matter properties must be described accurately at the same time,21,24,68 are well suited testing grounds for hybrid vdW-DFs. An additional benefit of hybrid vdW-DFs, not explored here, is the potential for accurate computations of excited-state properties. It is a practical and conceptual advantage to be able to compute the excited state properties of vdW bonded materials using the same hybrid vdW-DF version that was first used to relax the atomic coordinates.

The authors thank Zhenfei Liu and Elsebeth Schröder for insightful discussions. Work in Sweden is supported by the Swedish Research Council (VR), Chalmers Area-of-Advance Materials, and the Chalmers e-Science centre. K.B. acknowledges support from the Research Council of Norway (Nos. 250346 and 228854) and a NOTUR allocation. Work by J.H.L., T.R., and J.B.N. is supported by the Center for Gas Separations Relevant to Clean Energy Technologies, an Energy Frontier Research Center, funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0001015. Portions of this work were performed at the Molecular Foundry, supported by the Office of Science, Office of Basic Energy Sciences, U.S. Department of Energy, under Contract No. DE-AC02-05CH11231. Computational resources were provided by DOE (LBNL Lawrencium and NERSC) and by the Swedish National Infrastructure for Computing (SNIC), under Contract No. 2016-10-12, from the Chalmers Centre for computing, Science and Engineering (C3SE).

1.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
,
5648
(
1993
).
2.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
, “
Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields
,”
J. Phys. Chem.
98
,
11623
11627
(
1994
).
3.
C.
Adamo
and
V.
Barone
, “
Towards reliable density functional methods without adjustable parameters: The PBE0 model
,”
J. Chem. Phys.
110
,
6158
(
1999
).
4.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
, “
Hybrid functionals based on a screened Coulomb potential
,”
J. Chem. Phys.
118
,
8207
(
2003
).
5.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
, “
Erratum: ‘Hybrid functionals based on a screened Coulomb potential [J. Chem. Phys. 118, 8207 (2003)]
,’”
J. Chem. Phys.
124
,
219906
(
2006
).
6.
T.
Yanai
,
D. P.
Pew
, and
N. C.
Handy
, “
A new hybrid exchange-correlation functional using the Coulomb-attenuation method (CAM-B3LYP)
,”
Chem. Phys. Lett.
393
,
51
(
2004
).
7.
L.
Goerigk
and
S.
Grimme
, “
A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions
,”
Phys. Chem. Chem. Phys.
13
,
6670
6688
(
2011
).
8.
A. D.
Becke
, “
Perspective: Fifty years of density-functional theory in chemical physics
,”
J. Chem. Phys.
140
,
18A301
(
2014
).
9.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
789
(
1988
).
10.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
11.
K.
Burke
, “
Perspective on density functional theory
,”
J. Chem. Phys.
136
,
150901
(
2012
).
12.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
, “
Rationale for mixing exact exchange with density functional approximations
,”
J. Chem. Phys.
105
,
9982
9985
(
1996
).
13.
K.
Burke
,
M.
Ernzerhof
, and
J. P.
Perdew
, “
The adiabatic connection method: A non-empirical hybrid
,”
Chem. Phys. Lett.
265
,
115
120
(
1997
).
14.
M.
Ernzerhof
,
J. P.
Perdew
, and
K.
Burke
, “
Coupling-constant dependence of atomization energies
,”
Int. J. Quantum Chem.
64
,
285
295
(
1997
).
15.
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
, “
Van der Waals density functional for general geometries
,”
Phys. Rev. Lett.
92
,
246401
(
2004
).
16.
T.
Thonhauser
,
V. R.
Cooper
,
S.
Li
,
A.
Puzder
,
P.
Hyldgaard
, and
D. C.
Langreth
, “
Van der Waals density functional: Self-consistent potential and the nature of the van der Waals bond
,”
Phys. Rev. B.
76
,
125112
(
2007
).
17.
K.
Lee
,
È. D.
Murray
,
L.
Kong
,
B. I.
Lundqvist
, and
D. C.
Langreth
, “
Higher-accuracy van der Waals density functional
,”
Phys. Rev. B
82
,
081101
(
2010
).
18.
K.
Berland
and
P.
Hyldgaard
, “
Exchange functional that tests the robustness of the plasmon description of the van der Waals density functional
,”
Phys. Rev. B
89
,
035412
(
2014
).
19.
P.
Hyldgaard
,
K.
Berland
, and
E.
Schröder
, “
Interpretation of van der Waals density functionals
,”
Phys. Rev. B
90
,
075148
(
2014
).
20.
T.
Thonhauser
,
S.
Zuluaga
,
C. A.
Arter
,
K.
Berland
,
E.
Schröder
, and
P.
Hyldgaard
, “
Spin signature of nonlocal correlation binding in metal-organic frameworks
,”
Phys. Rev. Lett.
115
,
136402
(
2015
).
21.
K.
Berland
,
V. R.
Cooper
,
K.
Lee
,
E.
Schröder
,
T.
Thonhauser
,
P.
Hyldgaard
, and
B. I.
Lundqvist
, “
Van der Waals forces in density functional theory: A review of the vdW-DF method
,”
Rep. Prog. Phys.
78
,
066501
(
2015
).
22.
G.
Román-Pérez
and
J. M.
Soler
, “
Efficient implementation of a van der Waals density functional: Application to double-wall carbon nanotubes
,”
Phys. Rev. Lett.
103
,
096102
(
2009
).
23.
J.
Klimeš
and
A.
Michaelides
, “
Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory
,”
J. Chem. Phys.
137
,
120901
(
2012
).
24.
K.
Berland
,
C. A.
Arter
,
V. R.
Cooper
,
K.
Lee
,
B. I.
Lundqvist
,
E.
Schröder
,
T.
Thonhauser
, and
P.
Hyldgaard
, “
Van der Waals density functionals built upon the electron-gas tradition: Facing the challenge of competing interactions
,”
J. Chem. Phys.
140
,
18A539
(
2014
).
25.
S. D.
Chakarova-Käck
,
E.
Schröder
,
B. I.
Lundqvist
, and
D. C.
Langreth
, “
Application of van der Waals density functional to an extended system: Adsorption of benzene and naphthalene on graphite
,”
Phys. Rev. Lett.
96
,
146107
(
2006
).
26.
K.
Berland
,
S. D.
Chakarova-Käck
,
V. R.
Cooper
,
D. C.
Langreth
, and
E.
Schröder
, “
A van der Waals density functional study of adenine on graphene: Single-molecular adsorption and overlayer binding
,”
J. Phys.: Condens. Matter
23
,
135001
(
2011
).
27.
P.
Erhart
,
P.
Hyldgaard
, and
D.
Lindroth
, “
Microscopic origin of thermal conductivity reduction in disordered van der Waals solids
,”
Chem. Mater.
27
,
5511
(
2015
).
28.
W.
Liu
,
J.
Carrasco
,
B.
Santra
,
A.
Michaelides
,
M.
Scheffler
, and
A.
Tkatchenko
, “
Benzene adsorbed on metals: Concerted effect of covalency and van der Waals bonding
,”
Phys. Rev. B
86
,
245405
(
2012
).
29.
J.
Carrasco
,
B.
Santra
,
J.
Klimeš
, and
A.
Michaelides
, “
To wet or not to wet? Dispersion forces tip the balance for water ice on metals
,”
Phys. Rev. Lett.
106
,
026101
(
2011
).
30.
H.
Rydberg
,
M.
Dion
,
N.
Jacobson
,
E.
Schröder
,
P.
Hyldgaard
,
S. I.
Simak
,
D. C.
Langreth
, and
B. I.
Lundqvist
, “
Van der Waals density functional for layered structures
,”
Phys. Rev. Lett.
91
,
126402
(
2003
).
31.
E.
Londero
,
E. K.
Karlson
,
M.
Landahl
,
D.
Ostrovskii
,
J. D.
Rydberg
, and
E.
Schröder
, “
Desorption of n-alkanes from graphene: A van der Waals density functional study
,”
J. Phys.: Condens. Matter
24
,
424212
(
2012
).
32.
T.
Björkman
, “
Testing several recent van der Waals density functionals for layered structures
,”
J. Chem. Phys.
141
,
074708
(
2014
).
33.
H.
Peelaers
and
C. G.
Van de Walle
, “
First-principles study of van der Waals interactions in MoS2 and MoO3
,”
J. Phys.: Condens. Matter
26
,
305502
(
2014
).
34.
D. O.
Lindroth
and
P.
Erhart
, “
Thermal transport in van der Waals solids from first-principles calculations
,”
Phys. Rev. B
94
,
115205
(
2016
).
35.
C.
Ambrosch-Draxl
,
D.
Nabok
,
P.
Puschnig
, and
C.
Meisenbichler
, “
The role of polymorphism in organic thin films: Oligoacenes investigated from first principles
,”
New J. Phys.
11
,
125010
(
2009
).
36.
K.
Berland
and
P.
Hyldgaard
, “
Structure and binding in crystals of cagelike molecules: Hexamine and platonic hydrocarbons
,”
J. Chem. Phys.
132
,
134705
(
2010
).
37.
T.
Rangel
,
K.
Berland
,
S.
Sharifzadeh
,
F.
Brown-Altvater
,
K.
Lee
,
P.
Hyldgaard
,
L.
Kronik
, and
J. B.
Neaton
, “
Structural and excited-state properties of oligoacene crystals from first principles
,”
Phys. Rev. B
93
,
115206
(
2016
).
38.
J.
Klimeš
,
D. R.
Bowler
, and
A.
Michaelides
, “
Van der Waals density functionals applied to solids
,”
Phys. Rev. B
83
,
195131
(
2011
).
39.
L.
Gharaee
,
P.
Erhart
, and
P.
Hyldgaard
, “
Finite-temperature properties of non-magnetic transition metals: Comparison of the performance of constraint-based semi and nonlocal functionals
,”
Phys. Rev. B
95
,
085147
(
2017
).
40.
T.
Schwabe
and
S.
Grimme
, “
Double-hybrid density functionals with long-range dispersion corrections: Higher accuracy and extended applicability
,”
Phys. Chem. Chem. Phys.
9
,
3397
3406
(
2007
).
41.
J.-D.
Chai
and
M.
Head-Gordon
, “
Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections
,”
Phys. Chem. Chem. Phys.
10
,
6615
6620
(
2008
).
42.
S.
Grimme
, “
Density functional theory with London dispersion corrections
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
211
228
(
2011
).
43.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu
,”
J. Chem. Phys.
132
,
154104
(
2010
).
44.
A.
Tkatchenko
and
M.
Scheffler
, “
Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data
,”
Phys. Rev. Lett.
102
,
073005
(
2009
).
45.
A.
Tkatchenko
,
R. A.
DiStasio
,
R.
Car
, and
M.
Scheffler
, “
Accurate and efficient method for many-body van der Waals interactions
,”
Phys. Rev. Lett.
108
,
236402
(
2012
).
46.
N.
Marom
,
A.
Tkatchenko
,
M.
Rossi
,
V. V.
Gobre
,
O.
Hod
,
M.
Scheffler
, and
L.
Kronik
, “
Dispersion interactions with density-functional theory: Benchmarking semiempirical and interatomic pairwise corrected density functionals
,”
J. Chem. Theory Comput.
7
,
3944
(
2011
).
47.
O. A.
Vydrov
and
T.
Van Voorhis
, “
Nonlocal van der Waals density functional: The simpler the better
,”
J. Chem. Phys.
133
,
244103
(
2010
).
48.
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
, “
Erratum: ‘Van der Waals density functional for general geometries [Phys. Rev. Lett. 92, 246401 (2004)]
,’”
Phys. Rev. Lett.
95
,
109902
(
2005
).
49.
V. R.
Cooper
, “
Van der Waals density functional: An appropriate exchange functional
,”
Phys. Rev. B
81
,
161104
(
2010
).
50.
I.
Hamada
, “
Van der Waals density functional made accurate
,”
Phys. Rev. B
89
,
121103
(
2014
).
51.
O.
Gunnarsson
and
B. I.
Lundqvist
, “
Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism
,”
Phys. Rev. B
13
,
4274
4298
(
1976
).
52.
D. C.
Langreth
and
J. P.
Perdew
, “
Exchange-correlation energy of a metallic surface: Wave-vector analysis
,”
Phys. Rev. B
15
,
2884
2901
(
1977
).
53.
P.
Nozières
and
D.
Pines
, “
Correlation energy of a free electron gas
,”
Phys. Rev.
111
,
442
454
(
1958
).
54.
K.
Berland
and
P.
Hyldgaard
, “
Analysis of van der Waals density functional components: Binding and corrugation of benzene and C60 on boron nitride and graphene
,”
Phys. Rev. B
87
,
205421
(
2013
).
55.
J.
Schwinger
, “
Thomas-Fermi model: The second correction
,”
Phys. Rev. A
24
,
2353
2361
(
1981
).
56.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A. D.
Corso
,
S.
de Gironcoli
,
S.
Fabris
,
G.
Fratesi
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitsonen
,
A.
Smogunov
,
P.
Umari
, and
R. M.
Wentzcovitch
, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
57.
N.
Troullier
and
J. L.
Martins
, “
Efficient pseudopotentials for plane-wave calculations
,”
Phys. Rev. B
43
,
1993
2006
(
1991
).
58.
M.
Fuchs
and
M.
Scheffler
, “
Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory
,”
Comput. Phys. Commun.
119
,
67
98
(
1999
).
59.
P.
Jurecka
,
J.
Sponer
,
J.
Cerny
, and
P.
Hobza
, “
Benchmark database of accurate (MP2 and CCSD(T) complete basis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs
,”
Phys. Chem. Chem. Phys.
8
,
1985
1993
(
2006
).
60.
L.
Goerigk
and
S.
Grimme
, “
A general database for main group thermochemistry, kinetics, and noncovalent interactions—Assessment of common and reparameterized (meta-)GGA density functionals
,”
J. Chem. Theory Comput.
6
,
107
126
(
2010
).
61.
J. A.
Pople
,
M.
Head-Gordon
,
D. J.
Fox
,
K.
Raghavachari
, and
L. A.
Curtiss
, “
Gaussian–1 theory: A general procedure for prediction of molecular energies
,”
J. Chem. Phys.
90
,
5622
(
1989
).
62.

Results of the G1 set (subset of G2-1) for the non-hybrid spin vdW-DF were presented in the supplementary material of Ref. 20. These calculations were based on norm-conserving pseudopotentials rather than ultrasoft pseudopotentials and were compared to the values of the G161 set rather than the updated reference values of the G2 sets.60 

63.
W.
Hujo
and
S.
Grimme
, “
Performance of the van der Waals density functional VV10 and (hybrid)GGA variants for thermochemistry and noncovalent interactions
,”
J. Chem. Theory Comput.
7
,
3866
3871
(
2011
).
64.
A.
Hjorth Larsen
,
M.
Kuisma
,
J.
Löfgren
,
Y.
Pouillon
,
P.
Erhart
, and
P.
Hyldgaard
, “
libvdwxc: A library for exchange-correlation functionals in the vdW-DF family
,”
Modell. Simul. Mater. Sci. Eng.
25
,
065004
(
2017
); e-print arXiv:1703.06999 [cond-mat.mtrl-sci].
65.
J.
Řezáč
,
K. E.
Riley
, and
P.
Hobza
, “
Benchmark calculations of noncovalent interactions of halogenated molecules
,”
J. Chem. Theory Comput.
8
,
4285
4292
(
2012
).
66.
D. C.
Langreth
and
S. H.
Vosko
, “
Exact electron-gas response functions at high density
,”
Phys. Rev. Lett.
59
,
497
500
(
1987
).
67.
É. D.
Murray
,
K.
Lee
, and
D. C.
Langreth
, “
Investigation of exchange energy density functional accuracy for interacting molecules
,”
J. Chem. Theory Comput.
5
,
2754
2762
(
2009
).
68.
D. C.
Langreth
,
B. I.
Lundqvist
,
S. D.
Chakarova-Käck
,
V. R.
Cooper
,
M.
Dion
,
P.
Hyldgaard
,
A.
Kelkkanen
,
J.
Kleis
,
L.
Kong
,
S.
Li
,
P. G.
Moses
,
E.
Murray
,
A.
Puzder
,
H.
Rydberg
,
E.
Schröder
, and
T.
Thonhauser
, “
A density functional for sparse matter
,”
J. Phys.: Condens. Matter
21
,
084203
(
2009
).