# Semilocal and Hybrid Meta-Generalized Gradient Approximations Based on the Understanding of the Kinetic-Energy-Density Dependence

## Abstract

We present a global hybrid meta-generalized gradient approximation (meta-GGA) with three empirical parameters, as well as its underlying semilocal meta-GGA and a meta-GGA with only one empirical parameter. All of them are based on the new meta-GGA resulting from the understanding of kinetic-energy-density dependence [J. Chem. Phys. 137, 051101 (2012)]. The obtained functionals show robust performances on the considered molecular systems for the properties of heats of formation, barrier heights, and noncovalent interactions. The pair-wise additive dispersion corrections to the functionals are also presented.

## I Introduction

The Kohn-Sham (KS) density functional theory (DFT) (1); (2); (3) is one of the most widely used electronic structure theories for atoms, molecules and solids in various areas of physics, chemistry and molecular biology due to its computational efficiency and useful accuracy. It simplifies a many-electron wave-function problem to an auxiliary one-electron problem, with only its exchange-correlation part carrying the many-electron effects to be approximated in practice. Among numerous exchange-correlation approximations, the local spin density approximation (LSDA) (1); (4); (5), the standard Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) (6), and the Becke-3-Lee-Yang-Parr (B3LYP) (7); (8); (10); (9) hybrid GGA dominate the user market of DFT (11). The former two are efficient semilocal functionals widely used for extended systems, while the latter is a computationally more expensive nonlocal functional that hybridizes a GGA with the exact exchange energy and is popular for finite systems. At the semilocal level, however, meta-GGA (MGGA) is the highest rung of the so-called Jacob’s ladder in DFT (12) and potentially the most accurate one (13), which can also serve as a better base for hybridizing with the exact exchange energy.

Semilocal approximations (e.g., Refs. (4); (5); (6); (14); (15); (16)) of the form

(1) |

require only a single integral over real space and so are practical even for large molecules or unit cells. In Eq. (1), and are the electron densities of spin and , respectively, the local gradients of the spin densities, the kinetic energy densities of the occupied KS orbitals of spin (), and the approximate exchange-correlation energy per electron. All equations are in atomic units. In addition to and , the only ingredients in LSDA, GGAs also use the density gradients; MGGAs additionally include the kinetic energy densities . With the encoded information of shell structures from the kinetic energy densities, MGGAs can distinguish different orbital-overlap regions (17), and thus deliver simultaneous accurate ground-state properties for molecules, surfaces, and solids – which couldn’t be obtained with a GGA or LSDA (13). Computationally, MGGAs are not much more expensive than LSDA or GGA (18); (19). In computations for molecules containing transition-metal atoms, the Tao-Perdew-Staroverov-Scuseria (TPSS) (15) MGGA is only slower (19) than PBE.

To understand the performance of MGGAs, a recent study (17) investigated the effect of -dependence on MGGAs through the inhomogeneity parameter . Here, , is the von Weizscker kinetic energy density, and is the kinetic energy density of the uniform electron gas (UEG). is a lower bound on with only for a single-orbital region (20). Therefore, measures locally how far deviates from being of single-orbital character on the scale of the UEG, and thus characterizes the extent of orbital overlap. Besides the ability of to distinguish different orbital-overlap regions, Ref. (17) further showed that the effect of the -dependence on MGGAs is qualitatively equivalent to that of the dependence on the reduced density gradient , another dimensionless inhomogeneity parameter that measures how fast and how much the density varies on the scale of the local Fermi wavelength with . The -dependence is well understood at the GGA level (6); (14), which MGGAs inherit (15); (16); (17).

As a result of the study on the -dependence (17), a simple exchange functional, where the -dependence is disentangled from the -dependence, was constructed as an interpolation between the single-orbital regime () and the slowly-varying density regime () (17), which we will discuss more in Section III. In analogy to the monotonically increasing dependence that GGAs are often designed to have to favor more inhomogeneous electron densities, the simple exchange functional has a monotonically decreasing dependence to favor more the single-orbital regions. When combined with a variant of the PBE correlation (16); (17), the resulting MGGAs (denoted as MGGAMS with MS standing for ”made simple”) (17) performs equally well for atoms, molecules, surfaces, and solids, with an overall performance that is comparable to the sophisticated revised TPSS (revTPSS) MGGA (16). Moreover, MGGA_MS yields excellent binding energies for the W6 water clusters (17), even better than those from the highly parametrized M06-L MGGA (21) whose training sets include noncovalent interactions (hydrogen bonds and van der Waals interactions). Further tests on general main group thermochemistry, kinetic, and noncovalent interactions (22) showed that MGGA_MS, though not as good as M06-L, systematically improves the descriptions for the noncovalent interactions over conventional semilocal functionals, e.g., PBE, TPSS, and revTPSS.

By including training sets of noncovalent interactions, the M06-L MGGA was trained to capture medium-range exchange and correlation energies that dominate equilibrium structures of noncovalent complexes (21). For example, M06-L delivers good binding energies for the S22 set that is dominated by noncovalent interactions and will be discussed more in Sec. IV, good structural and energetical properties of layered materials bound by van der Waals interactions (e.g., graphite, hexagonal boron nitride, and molybdenum disulfite (23)), and good adsorption energies of aromatic molecules adsorbed on coinic metal surfaces (24). M06-L represents one of the developments of MGGAs that are fitted to a variety of training sets with numerous parameters to try to improve different properties together (21); (25); (26). However, it should be stressed that too many fitting parameters could cause problems, e.g., convergence of calculations and the related smoothness of binding curves (27), and the prediction for systems far from training sets might not be realistic (28).

Semilocal (sl) approximations can be reasonably accurate for the near-equilibrium and compressed ground-state properties of ââordinaryââ matter, where neither strong correlation nor long-range van der Waals interaction is important. However, stretched bonds that arise in transition states of chemical reactions or in the dissociation limits of some radical or heteronuclear molecules, etc., require full nonlocality (29). The global hybrid functionals discussed later can provide such needed nonlocality. Long-range van der Waals interactions, originated from spontaneous correlations between two distant electron densities, also require full nonlocality (30), which is different from the previous one (31). To obtain long-range van der Waals interactions, there are various DFT-based dispersion techniques, including the post-DFT empirical pairwise potential corrections (32); (33); (34); (35); (36); (37); (38); (39) and the van der Waals density functionals (vdW-DFs) (30); (40); (41), which are based on the pairwise additivity and problematic for metallic systems (42); (43). See Ref. (44) for an overview of the progess on DFT-based dispersion techniques. Conventional semilocal functionals, e.g., the widely used PBE GGA, even lack the ability to describe noncovalent interactions near equilibrium, while the M06-L MGGA and MGGA_MS as well as its variant proposed here improve over PBE.

The global hybrid (gh) idea, due to Becke (7); (45), introduces some full nonlocality into the calculation but only at the level of , which can be evaluated semianalytically from the Kohn-Sham orbitals in some computer codes. In its simplest (one parameter) version (46):

(2) |

where the exact-exchange mixing parameter takes an empirical value in the range . A rough but convincing argument (47) for why Eq. (2) works is: semilocal exchange-correlation typically overestimates atomization energies and underestimates energy barriers of chemical reactions, while exact exchange without correlation makes errors of the opposite sign, so that mixing of the two will yield better atomization energies and barriers than either alone. Therefore, values of vary with the choice of semilocal functionals when fitting Eq. (2) to formation enthalpies. The more a semilocal functional overbinds molecules, the larger the value of is. Csonka et al. (47) found to be 0.60 for the PBEsol GGA, 0.32 for the PBE GGA, and 0.1 for the TPSS and revTPSS MGGAs. Improving the underlying semilocal functional reduces the value of needed to fit the atomization energies or enthalpies of formation, which however can worsen the barrier heights. Some global hybrids (48); (49) thus employ highly fitted semilocal parts intended not to be accurate by themselves but to work well with a fraction of exact exchange.

The popular B3LYP (7); (8); (10); (9) hybrid GGA has a somewhat more complicated form and reads:

(3) | |||||

where (50) is the LSDA correlation energy in the random phase approximation which does not yield the correct uniform electron gas limit. In the equation, a=0.2, b=0.72, and c=0.81. If we use the similar idea as in Eq. (2) to define the underlying semilocal functional for B3LYP by replacing with , we end up with a GGA, denoted as BLYP, . BLYP differs from BLYP in reducing the gradient correction of the B88 exchange to (a+b)=92% and mixing the VWN3 and LYP correlations. This difference confirms the rough argument of the previous paragraph. Typically, BLYP underbinds molecules, and inclusion of exact exchange worsens the underbinding. The reduction of the gradient correction of the B88 exchange and the mixture of the VWN3 and LYP correlations then turn BLYP into an overbinding functional suitable for hybridization with the exact exchange. The idea of tuning the gradient correction of B3LYP can be transfered to the construction of our hybrid MGGA, where both the dependence and the dependence are tuned.

## Ii Computational details

The functionals described here were implemented in the development version of the Gaussian electronic structure program (51). All calculations employ the fully uncontracted 6-311++G() basis set (52); (53) to obtain benchmark quality results. We used the UltraFine grid with 99 radial shells and 590 angular points for numerical integration of the DFT XC potential. The Gaussian program was used for all calculations presented in this study.

Different training sets were used to adjust empirical parameters. For the AE6 and BH6 tests sets (54) we employed reference data from highly accurate quantum-chemical calculations (55). Experimental reference values were used for G2/97 (148 molecules) (56); (57) and BH42/03 (21 forward and reverse hydrogen transfer barrier heights) (58). The fitted parameters would not change much if the highly accurate quantum-chemical reference data from Ref. (59) would have been used.

The superset of G2/97 and G3-3 (75 molecules) (60) (a total of 223 molecules comprising the G3/99 test set) (61), the S22 test set for weak interactions (62), as well as the barrier height test sets HTBH38/04 (19 forward and reverse hydrogen transfer barrier heights) (58) and NHTBH38/04 (19 forward and reverse non-hydrogen transfer barrier heights) (63) are employed for assessment of the fitted functionals. The employed geometries and reference values are available from the Supporting Information of Refs. (64) and (62). The counterpoise correction (65) was used to reduce the basis set superposition errors for calculations of weak interactions. A vibrational scale factor of 0.9854 was used for calculations of heats of formation on the B3LYP/6-31G() level of theory as recommended in Ref. (66).

The performance of our fitted functionals will be compared to related and popular density functionals: PBE (6), TPSS (15), revTPSS (16), M06-L (21), MGGA_MS0 (17), PBEh (46), TPSSh (67), revTPSSh (16), B3LYP (8); (10); (9), BLYP (8); (10), and M06 (68). Calculations of open-shell species were carried out via spin-unrestricted formalisms. Errors are reported as .

## Iii Construction of semilocal and hybrid meta-GGAs

The semilocal exchange energy of a spin-unpolarized density can be written as:

(4) |

Here, is the exchange energy per particle of a UEG with , , and the enhancement factor with for LSDA. The expression for the spin-polarized case follows from the exact spin scaling (69). In Ref. (17), a simple exchange enhancement factor that disentangles the - and -dependence was proposed,

(5) |

where and . interpolates between and through . Here, in order to tune the dependence, we add a parameter in the denominator of the original , which reads now . For a slowly varying density where , is controlled by its numerator and we recover the second order gradient expansion with the first principle coefficient (70) as in Ref. (17). The parameter is determined for each to reproduce the exchange energy of the hydrogen atom, where . Since is negligible for and the numerator of modulates the behavior of for , mainly controls the behavior for . On the other hand, mainly controls the behavior of for large . In the original form (17), , and was determined by the Hartree-Fock exchange energy of the 12-noninteracting-electron hydrogenic anion with nuclear charge Z=1 (71). Fig. 2 of Ref. (17) showed that the shell regions of the hydrogenic anion are associated with and the intershell regions with , which thus can provide information to guide the exchange functional from to . Here, we release the constraint of the hydrogenic anion and rely on some molecular training sets to provide information to determine both and by fitting them to heats of formation and barrier heights. For the correlation part, we use the variant of the PBE correlation (16); (17).

When we plug the above exchange and correlation parts into Eq. (2), we obtain a global hybrid MGGA with three parameters, , , and , which control the dependence of the global hybrid on the density gradient, the kinetic energy density, and the exact exchange, respectively. In the next section, we describe the parametrization procedure and the performance of our best functionals.

## Iv Results and Discussions

To get a broad overview about the errors which the semilocal functional class MGGA_MS makes using different values of and , we define the average mean absolute error (AMAE) using the mean absolute error (MAE) values from the AE6 and BH6 test sets:

(6) |

and present the plotted surface of AMAE values in figure 1. The figure shows that there are two useful possibilities of fitting a semilocal functional of our MGGA_MS form: First, we can design a one-parameter functional by setting the parameter to one and fitting to minimize a chosen error measure. Second, both parameters can be fitted independently. The second choice will yield a lower AMAE value. For the independent fit, we also include the exact exchange admixture parameter into the fitting procedure. The corresponding semilocal functional will be obtained by setting to zero while keeping the other global hybrid parameters fixed.

Our chosen error measure contains the mean absolute errors (MAE’s) of the heats of formation of the G2/97 data set and barrier heights of the BH42/03 test sets:

(7) |

The resulting parameter sets, ME, and MAE values are shown in table 1 and compared to other popular semilocal and hybrid functionals. MGGA_MS0 is the original MGGA_MS functional from Ref. (17) with no parameter fitting to training sets while MGGA_MS1 is our semilocal one-parameter functional. MGGA_MS2h is the new hybrid functional and MGGA_MS2 its semilocal complement. Our hybrid MGGA_MS2h has a rather low exact exchange admixture of 9% as expected when a good semilocal functional is hybridized (47). Comparing the parameter set of MGGA_MS2 with the AMAE surface in figure 1 shows that its parameter set is quite close to the global minimum of the AMAE, although it has not been refitted after exact exchange admixture was reduced from 9% to 0%. The functional BLYP* is not commonly used, but it is the underlying semilocal functional of B3LYP. It is included in our tables for comparison purposes only. As shown in tables 1 and 2, BLYP* turns the underestimations of BLYP for atomization energies into overestimations. From all semilocal functionals in table 1, MGGA_MS1 performs best for G2/97 while M06-L performs best for BH42/03. MGGA_MS0, MGGA_MS1, and MGGA_MS2 functionals are the next best performers for the BH42/03 test set, which behave similarly despite very different parameter sets. For the hybrids, B3LYP performs best for G2/97 and M06 for BH42/03. Our global hybrid MGGA_MS2h performs similarly to B3LYP and PBEh for the BH42/03 test set and similar to PBEh for the G2/97 test set. Note that PBE, TPSS, revTPSS, and MGGA_MS0 are nonempirical functionals, while the number of empirical parameters – as indicated in the parentheses – of MGGA_MS1 (1), MGGA_MS2 (2), TPSSh (1), revTPSSh (1), and MGGA_MS2h (3) are an order of magnitude smaller than those of M06-L and M06.

MGGA_MS1 performs better than MGGA_MS2 and MGGA_MS0 for all molecular test sets discussed so far. This and figure 1 indicate that one can construct another MGGA_MS functional which performs even better than MGGA_MS1. However, figure 1 suggests that the improvement over MGGA_MS1 will not be significant. We found that , , and is close to optimal for the error measure , but the minimum in parameter space is rather shallow and broad. With MAE values for G2/97 of 4.5 kcal/mol and for BH42/03 of 5.6 kcal/mol, this variant would improve over MGGA_MS1 only slightly for G2/97 and not at all for BH42/03 on the expense of an additional empirical parameter. Thus, we drop this variant from now on.

G2/97 | BH42/03 | |||||||

ME | MAE | ME | MAE | |||||

PBE | 0.00 | 0.804 | – | – | -16.1 | 16.9 | -9.7 | 9.7 |

BLYP | 0.00 | – | – | – | -0.6 | 7.4 | -8.0 | 8.0 |

BLYP* | 0.00 | – | – | – | -8.1 | 9.0 | -8.4 | 8.4 |

TPSS | 0.00 | 0.804 | – | – | -5.7 | 6.4 | -8.4 | 8.4 |

revTPSS | 0.00 | 0.804 | – | – | -4.7 | 5.7 | -7.7 | 7.7 |

M06-L | 0.00 | – | – | – | -3.0 | 5.1 | -4.4 | 4.4 |

MGGA_MS0 | 0.00 | 0.29 | 0.28771 | 1.0 | -0.6 | 6.8 | -5.9 | 6.0 |

MGGA_MS1 | 0.00 | 0.404 | 0.18150 | 1.0 | 2.3 | 4.9 | -5.5 | 5.6 |

MGGA_MS2 | 0.00 | 0.504 | 0.14601 | 4.0 | -0.8 | 5.1 | -5.8 | 5.8 |

PBEh | 0.25 | 0.804 | – | – | -2.6 | 5.0 | -4.7 | 4.7 |

TPSSh | 0.10 | 0.804 | – | – | -1.9 | 4.4 | -6.6 | 6.6 |

revTPSSh | 0.10 | 0.804 | – | – | -1.2 | 4.5 | -6.0 | 6.0 |

B3LYP | 0.20 | – | – | – | 0.9 | 3.1 | -4.7 | 4.7 |

M06 | 0.27 | – | – | – | -3.4 | 4.1 | -2.3 | 2.3 |

MGGA_MS2h | 0.09 | 0.504 | 0.14601 | 4.0 | 2.7 | 4.9 | -4.5 | 4.6 |

G3-3 | G3/99 | HTBH38/04 | NHTBH38/04 | |||||

ME | MAE | ME | MAE | ME | MAE | ME | MAE | |

PBE | -32.6 | 32.6 | -21.6 | 22.2 | -9.7 | 9.7 | -8.5 | 8.6 |

BLYP | 12.5 | 14.1 | 3.8 | 9.7 | -7.9 | 7.9 | -8.7 | 8.8 |

BLYP* | -5.7 | 8.6 | -7.3 | 8.8 | -8.4 | 8.4 | -8.7 | 8.8 |

TPSS | -6.4 | 6.5 | -6.0 | 6.5 | -8.2 | 8.2 | -9.0 | 9.1 |

revTPSS | -4.2 | 4.7 | -4.5 | 5.4 | -7.4 | 7.4 | -9.3 | 9.3 |

M06-L | -4.2 | 7.4 | -3.4 | 5.8 | -4.5 | 4.5 | -3.2 | 3.7 |

MGGA_MS0 | -5.9 | 11.8 | -2.4 | 8.5 | -5.6 | 5.7 | -6.4 | 6.9 |

MGGA_MS1 | 3.2 | 5.2 | 2.6 | 5.0 | -5.4 | 5.5 | -6.3 | 6.7 |

MGGA_MS2 | -4.1 | 7.6 | -1.9 | 6.0 | -5.7 | 5.7 | -7.0 | 7.3 |

PBEh | -9.6 | 10.4 | -4.9 | 6.8 | -4.7 | 4.7 | -3.2 | 3.7 |

TPSSh | -1.0 | 3.5 | -1.6 | 4.1 | -6.4 | 6.4 | -6.9 | 7.0 |

revTPSSh | 0.6 | 3.6 | -0.6 | 4.2 | -5.8 | 5.8 | -7.1 | 7.2 |

B3LYP | 7.9 | 8.2 | 3.3 | 4.8 | -4.5 | 4.6 | -4.6 | 4.7 |

M06 | -4.4 | 5.4 | -3.7 | 4.5 | -2.3 | 2.3 | -1.7 | 2.2 |

MGGA_MS2h | 1.5 | 5.1 | 2.3 | 4.9 | -4.4 | 4.4 | -5.3 | 5.7 |

To test the deviation from the exact exchange energy of the 12-noninteracting-electron hydrogenic anion, we calculate the exchange energies for the new parameter sets of (b, ), which are -1.8562 Ha for (b=1, ) of MGGA_MS1 and -1.8558 Ha for (b=4, ) of MGGA_MS2, respectively. In comparison with the exact one, -1.8596 Ha, the deviations are small, which indicates that the information of the hydrogenic anion is preserved by Eq. 5 and complemented by that from training sets. Figure 2 shows the exchange enhancement factors for different and values of our new functionals and compares them to the revTPSS functional. Figure 2 shows revTPSS has the order of limits anomaly at small and small (15); (16), which can be eliminated by regularizing revTPSS at these regions (72). The family of MGGA_MS doesn’t have such a problem by construction. The enhancement factors of different MGGAs considered here are on top of each other around at the curves of , due to the constraint of the second order gradient expansion. Figure 2 shows that each MGGA is insensitive to for large . For the curves, which is modulated by the constraint of the exchange energy of the hydrogen atom, gets smaller for small regions when a larger is used. The behavior of these MGGAs at large s is determined by the value of . The larger , the larger enhancement factor at large . Due to the qualitative equivalence between the and dependence (17), we observe for different MGGAs that the faster the increase of F with , the slower the decrease of F with .

It is not too surprising that our functionals fitted to G2/97 and BH42/03 are among the best performers for these test sets. Table 2 shows the performance for the heats of formation of the test sets G3-3 and G3/99 as well as the barrier heights of the data sets HTBH38/04 and NHTBH38/04. The Minnesota functionals perform best for the barrier heights, which is not too surprising as both barrier height test sets are contained in the training set of M06 and M06-L. The next best semilocal functional for both barrier height test sets is MGGA_MS1. Among the global hybrids, PBEh and B3LYP perform better than MGGA_MS2h for the NHTBH38/04 test set but not for the HTBH38/04 test set which was part of the training set. TPSSh performs best for G3-3 and G3/99 within our choice of global hybrids. Interestingly, all global hybrids in table 2 (except PBEh) perform very similar on average for the G3/99 test set. Among the semilocal functionals revTPSS performs best for G3-3 and MGGA_MS1 for G3/99.

Aside from heats of formation and barrier heights, we would like to test our functionals for a property which was not part of the training set. We choose the weak interactions of the S22 test set. Table 3 presents the ME and MAE values for the S22 test set and its subsets WI8 (including 8 dispersion-bound complexes), MI7 (including 7 mixed complexes), and HB7 (including 7 hydrogen-bonded complexes). Before discussing the results, we would like to test the convergence of MGGAs with respect to the integration grid for dispersion-bound complexes. Johnson et al (73) found that very large integration grids (much finer than the UltraFine) are required to remove oscillations in potential energy surfaces (PES) for dispersion-bound complexes for most of the MGGAs they studied. For the chosen simple dispersion-bound NeAr, we confirmed that the M06L MGGA, which was in Johnson’s test list, needs a very large grid to yield a smooth PES. However, our MGGAs converge much faster for the total energies than M06L does, and yield converged smooth curves with the UltraFine grid. The Minnesota functionals, which have been fitted for weak interactions, perform best for all these data sets of weak interactions. Among semilocal functionals, MGGA_MS0 is the second best performer for the S22 set and its subsets. MGGA_MS2 on average is comparable to MGGA_MS0 for the WI8 and MI7 subsets, but worse for the HB7 subset, while MGGA_MS1 loses the accuracy for the weak interactions compared to the other two in the MGGA_MS family. As expected for WI8 and MI7, where dispersion interactions are important, PBE, TPSS, and revTPSS yield large MAEs, which are worse than the other semilocal functionals considered here.

WI8 | MI7 | HB7 | S22 | |||||

ME | MAE | ME | MAE | ME | MAE | ME | MAE | |

PBE | 4.8 | 4.8 | 2.0 | 2.0 | 1.1 | 1.1 | 2.8 | 2.8 |

BLYP | 7.7 | 7.7 | 3.8 | 3.8 | 3.2 | 3.2 | 5.0 | 5.0 |

BLYP* | 7.3 | 7.3 | 3.6 | 3.6 | 2.8 | 2.8 | 4.7 | 4.7 |

TPSS | 6.0 | 6.0 | 2.7 | 2.7 | 1.9 | 1.9 | 3.7 | 3.7 |

revTPSS | 4.3 | 4.3 | 2.4 | 2.4 | 2.0 | 2.0 | 2.9 | 2.9 |

M06-L | 1.0 | 1.0 | 1.0 | 1.0 | 0.7 | 0.7 | 0.9 | 0.9 |

MGGA_MS0 | 3.0 | 3.0 | 1.2 | 1.2 | 0.8 | 0.8 | 1.7 | 1.7 |

MGGA_MS1 | 3.9 | 3.9 | 1.7 | 1.7 | 1.8 | 1.8 | 2.5 | 2.5 |

MGGA_MS2 | 3.0 | 3.0 | 1.3 | 1.3 | 1.4 | 1.4 | 1.9 | 1.9 |

PBEh | 4.6 | 4.6 | 1.8 | 1.8 | 0.8 | 0.9 | 2.5 | 2.5 |

TPSSh | 5.8 | 5.8 | 2.6 | 2.6 | 1.8 | 1.8 | 3.5 | 3.5 |

revTPSSh | 5.0 | 5.0 | 2.3 | 2.3 | 1.8 | 1.8 | 3.1 | 3.1 |

B3LYP | 6.5 | 6.5 | 3.0 | 3.0 | 2.0 | 2.0 | 4.0 | 4.0 |

M06 | 1.4 | 1.4 | 0.9 | 0.9 | 0.7 | 0.7 | 1.0 | 1.0 |

MGGA_MS2h | 3.0 | 3.0 | 1.3 | 1.3 | 1.2 | 1.2 | 1.9 | 1.9 |

MGGA_MS0+D3 | MGGA_MS1+D3 | MGGA_MS2+D3 | MGGA_MS2h+D3 | |||||

test set | ME | MAE | ME | MAE | ME | MAE | ME | MAE |

G2/97 | -0.9 | 7.0 | 1.8 | 5.0 | -1.1 | 5.3 | 2.3 | 4.9 |

BH42/03 | -6.1 | 6.1 | -5.9 | 5.9 | -6.0 | 6.0 | -4.7 | 4.7 |

G3-3 | -7.4 | 13.1 | 0.6 | 5.7 | -5.6 | 9.0 | -0.0 | 5.8 |

G3/99 | -3.1 | 9.1 | 1.4 | 5.3 | -2.6 | 6.6 | 1.6 | 5.2 |

HTBH38/04 | -5.8 | 5.9 | -5.8 | 5.8 | -6.0 | 6.0 | -4.6 | 4.6 |

NHTBH38/04 | -6.5 | 7.0 | -6.5 | 6.9 | -7.1 | 7.4 | -5.4 | 5.8 |

WI8 | -0.0 | 0.4 | -0.2 | 0.4 | -0.2 | 0.3 | -0.1 | 0.3 |

MI7 | -0.3 | 0.3 | -0.3 | 0.3 | -0.2 | 0.3 | -0.2 | 0.3 |

HB7 | -0.3 | 0.4 | 0.2 | 0.2 | 0.2 | 0.3 | 0.1 | 0.2 |

S22 | -0.2 | 0.3 | -0.1 | 0.3 | -0.0 | 0.3 | -0.1 | 0.2 |

Ref. (44) presents a collection of different techniques for functionals sorted by the way weak interactions are accounted for: step 1 (DFT-D (37), etc), step 2 (DFT-D3 (38), vdW(TS) (39), etc), step 3 (long-range density functionals (30); (40); (41)), and step 4 and higher (RPA (74), etc). Representative functionals of each group are tested for the S22 test set in Ref. (44). Note that the basis sets differ compared to our calculations, but the results can be compared roughly. MGGA_MS0 and MGGA_MS2 yield results that are close to the vdW-DF (30) which is on step 3 of this classification, while MGGA_MS1 fits best to the ground (below step 1). It should be noted that the classification is methodological and not based on performance as some functionals in step 1 and 2 outperform some funtionals on step 3 in Ref. (44). Table 4 gives results of the MGGA_MS functionals with the D3 correction of Grimme et al (38). The parameters in the D3 correction term were optimized by fitting to the S22 data set for each functional. We found that the optimization is insensitive with respect to the parameter , and thus set . Then, there is only one fitting parameter left in the D3 correction term. , which is the scaling factor of cutoff radii and determines the range of the dispersion correction inversely (38), is 1.15, 1.05, 1.14, and 1.14 for MGGA_MS0, MGGA_MS1, MGGA_MS2, and MGGA_MS2h, respectively. All the four MGGA_MS functionals with the D3 corrections deliver for the S22 data set similar MAEs around 0.3 kcal/mol, which are among the best of the D3-corrected functionals tested in Ref. (38). The smaller for MGGA_MS1 indicates that MGGA_MS1 captures less intermediate-range weak interactions than the other three do, as also shown in Table 3. Table 4 also shows that the D3 corrections have noticeable effects on the formation energies and the barrier heights. MGGA_MS1+D3 is the best among the three MGGAs, while the hybrid MGGA_MS2h+D3 improves further the formation energies and the barrier heights.

The functional pairs of global hybrid and its corresponding semilocal complement perform very similarly for weak interactions on average. The pair of BLYP* and B3LYP is the outlier, which might be related to the fact that B88 is too repulsive for rare gas dimers even when comparing to the exact exchange (75). This general trend suggests that the performance of a hybrid of a semilocal functional inherits its performance for weak interactions from its corresponding pure semilocal functional. A possible explanation is that weak interactions are mainly correlation effects where exact exchange does not contribute. On the other hand, heats of formation, atomization energies, and barrier heights are significantly affected by correlation and exact exchange contributions. However, in our constructions, MGGA_MS2 is a byproduct of MGGA_MS2h, which is fitted to the atomization energies of the G2/97 test set and the barrier heights of the BH42/03. Surprisingly, MGGA_MS2h performs reasonably well for the S22 set, much better than the other four hybrids, namely, PBEh, TPSSh, revTPSSh, and B3LYP. Note the improvement of MGGA_MS2h over the other four hybrids on the weak interactions is not a result of sacrifice of accuracy for atomization energies and barrier heights as shown in Table 2. The good performance of MGGA_MS2h on the S22 set is then transfered to MGGA_MS2.

When considering both Tables 2 and 3 for the properties of atomization energies, barrier heights, and weak interactions, our three-parameter global hybrid MGGA_MS2h overall is the second best hybrid after M06 that is heavily parameterized. Compared to MGGA_MS2h, our semilocal functional MGGA_MS1 performs surprisingly well for atomization energies of molecules and barrier heights of chemical reactions although it contains just one empirical parameter. However, the gain is obtained at the price of descriptions for weak interactions as demonstrated by the comparison with MGGA_MS0 in Table 3. This shows the limit of tuning only the dependence of the enhancement factor. Fortunately, by tuning the dependence additionally, MGGA_MS2 improves over MGGA_MS0 for the atomization energies significantly while remaining comparable for the barrier heights and weak interactions.

## V Conclusions

By taking advantage of the understanding on the kinetic-energy-density dependence of MGGAs (17), we construct a global hybrid meta-generalized gradient approximation (hybrid meta-GGA) with three empirical parameters (MGGA_MS2h), which has robust performances on the molecular systems for the properties of heats of formation, barrier heights, and noncovalent interactions. The derived underlying MGGA_MS2 improves over the original MGGA_MS0 (17) for heats of formation significantly while retaining good performance for barrier heights and weak interactions. By relaxing the parameter of MGGA_MS0, we also obtained the one-parameter functional MGGA_MS1 that performs surprisingly well for heats of formations and barrier heights.

###### Acknowledgements.

Most calculations were done at Rice University. This work was supported by the National Science Foundation under CHE-1110884 and the Welch Foundation (C-0036). RH thanks the Deutsche Forschungsgemeinschaft (DFG grant HA 5711/2-1). JS, BX, and JPP acknowledge NSF support under Grant No. DMR-0854769, and Cooperative Agreement No. EPS-1003897, with additional support from the Louisiana Board of Regents.### References

- W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).
- R.G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford University Press, Oxford, 1989).
- J.P. Perdew and S. Kurth, in A Primer in Density Functional Theory (Springer Lecture Notes in Physics, Vol. 620, 2003).
- J.P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
- J. Sun, J.P. Perdew and M. Seidl, Phys. Rev. B 81, 085123 (2010).
- J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
- A.D. Becke, Phys. Rev. A 38, 3098 (1988).
- P.J. Stephens, F.J. Devlin, C.F. Chabalowski, and M.J. Frisch, J. Phys. Chem. 98, 11623 (1994).
- C. Lee, W. Yang, and R.G. Parr, Phys. Rev. B 37, 785 (1988).
- K. Burke, J. Chem. Phys. 136, 150901 (2012).
- J. P. Perdew and K. Schmidt, in Density Functional Theory and Its Applications to Materials, edited by V. E. van Doren, C. van Alsenoy, and P. Geerlings (American Institute of Physics, 2001).
- J. Sun, M. Marsman, A. Ruzsinszky, G. Kresse, and J.P. Perdew, Phys. Rev. B 83, 121410(R) (2011).
- J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
- J. Tao, J.P. Perdew, V.N. Staroverov, and G.E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
- J.P. Perdew, A. Ruzsinszky, G.I. Csonka, L.A. Constantin, and J. Sun, Phys. Rev. Lett. 103, 026403 (2009); ibid. 106, 179902 (2011) (E).
- J. Sun, B. Xiao, and A. Ruzsinszky, J. Chem. Phys. 137, 051101 (2012).
- V.N. Staroverov, G.E. Scuseria, J. Tao, and J.P. Perdew, J. Chem. Phys. 119, 12129 (2003).
- F. Furche and J.P. Perdew, J. Chem. Phys. 124, 044103 (2006).
- S. Kurth, J.P. Perdew, and P. Blaha, Int. J. Quantum Chem. 75, 889 (1999).
- Y. Zhao and D.G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
- P. Hao, J. Sun, B. Xiao, A. Ruzsinszky, G.I. Csonka, J. Tao, and J.P. Perdew, J. Chem. Theor. Comput. 10, 1021 (2012).
- G.K.H. Madsen, L. Ferrighi, and B. Hammer, J. Phys. Chem. Lett. 1, 515 (2010).
- L. Ferrighi, G.K.H. Madsen, and B. Hammer, J. Chem. Phys. 135, 084704 (2011).
- R. Peverati and D.G. Truhlar, J. Phys. Chem. Lett. 3, 117 (2012).
- T. Van Voorhis and G. E. Scuseria, J. Chem. Phys., 109, 400 (1998).
- L. Goerigk and S. Grimme, Phys. Chem. Chem. Phys. 13, 6670 (2011).
- J.C. Snyder, M. Rupp, K. Hansen, K-R. Mueller, and K. Burke, Phys. Rev. Lett. 108, 253002 (2012).
- J.P. Perdew, V.N. Staroverov, J. Tao, and G.E. Scuseria, Phys. Rev. A 78, 052513 (2008).
- M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, and B.I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
- A. Ruzsinszky, J.P. Perdew, and G.I. Csonka, J. Chem. Theor. Comput. 6, 127 (2010).
- Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002).
- A.D. Becke, J. Chem. Phys. 122, 064101 (2005).
- A.D. Becke, J. Chem. Phys. 122, 154104 (2005).
- A.D. Becke, J. Chem. Phys. 123, 024101 (2005).
- A.D. Becke, J. Chem. Phys. 123, 154101 (2005).
- S. Grimme, J. Comput. Chem. 27, 1787 (2006).
- S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
- A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
- O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. 103, 063004 (2009).
- J. Klimeš, D. R. Bowler, and A. Michaelides, J. Phys.: Condens. Matter 22, 022201 (2010).
- J.F. Dobson, A. White, and A. Rubio, Phys. Rev. Lett. 96, 073201 (2006).
- A. Ruzsinszky, J.P. Perdew, J. Tao, G.I. Csonka, and J.M. Pitarke, Phys. Rev. Lett. (to appear).
- J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
- A.D. Becke, J. Chem. Phys. 98, 1372 (1993).
- J.P. Perdew, K. Burke, and M. Ernzerhof, J. Chem. Phys. 105, 9982 (1996).
- G.I. Csonka, J.P. Perdew, and A. Ruzsinszky, J. Chem. Theor. Comput. 6, 3688 (2010).
- A.D. Becke, J. Chem. Phys. 107, 8554 (1997).
- Y. Zhao and D.G. Truhlar, Acc. Chem. Res. 41, 157 (2008).
- S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
- Frisch, Trucks, Schlegel, et al., Gaussian development version revision g.01, Gaussian Inc., Wallingford, CT (2009).
- A.D. McLean and G.S. Chandler, J. Chem. Phys. 72, 5639 (1980).
- R. Krishnan, J.S. Binkley, R. Seeger, and J.A. Pople, J. Chem. Phys. 72, 650 (1980).
- B.J. Lynch, D.G. Truhlar, J. Phys. Chem. A 107, 8996 (2003); erratum: 108, 1460 (2004).
- R. Haunschild and W. Klopper, Theor. Chem. Acc. 131, 1112 (2012).
- L.A. Curtiss, K. Raghavachari, P.C. Redfern, J.A. Pople, J. Chem. Phys. 106, 1063 (1997).
- L.A. Curtiss, P.C. Redfern, K. Raghavachari, J.A. Pople, J. Chem. Phys. 109, 42 (1998).
- Y. Zhao, B.J. Lynch, and D. G. Truhlar, J. Phys. Chem. A 108, 2715 (2004).
- R. Haunschild and W. Klopper, J. Chem. Phys. 136, 164102 (2012).
- L.A. Curtiss, K. Raghavachari, P.C. Redfern, J.A. Pople, J. Chem. Phys. 112, 7374 (2000).
- COF was included in the G2/97 and G3/99 test sets following Ref. (67).
- P. Jurečka, J. Šponer, J. Černý, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006).
- Y. Zhao, N. Gónzales-García, and D.G. Truhlar, J. Phys. Chem. A 109, 2012 (2005); erratum: 109, 4942 (2006).
- R. Haunschild, B.G. Janesko, and G.E. Scuseria, J. Chem. Phys. 131, 154112 (2009).
- S.F. Boys and F. Bernadi, Mol. Phys. 19, 553 (1970).
- L.A. Curtiss, P.C. Redfern, K. Raghavachari, and J.A. Pople, J. Chem. Phys. 114, 108 (2001).
- V.N. Staroverov, G.E. Scuseria, J. Tao, and J.P. Perdew, J. Chem. Phys. 119, 12129 (2003); erratum: 121, 11508 (2004).
- Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008).
- J.P. Perdew, S. Kurth, A. Zupan, and P. Blaha, Phys. Rev. Lett. 82, 2544 (1999); 82, 5179(E) (1999).
- P.R. Antoniewicz and L. Kleinman, Phys. Rev. B 31, 6779 (1985).
- V.N. Staroverov, G.E. Scuseria, J.P. Perdew, J.Tao, and E.R. Davidson, Phys. Rev. A 70, 012502 (2004)
- A. Ruzsinszky, J. Sun, B. Xiao, and G.I. Csonka, J. Chem. Theo. Comp. 8, 2078 (2012).
- E.R. Johnson, A.D. Becke, C.D. Sherrill, and G.A. DiLabio, J. Chem. Phys. 131, 034111 (2009).
- F. Furche, Phys. Rev. B 64, 195120 (2001).
- F.O. Kannemann and A.D. Becke, J. Chem. Theor. Comput. 5, 721 (2009).