New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Doubly hybrid density functional for accurate descriptions of nonbond interactions, thermochemistry, and thermochemical kinetics

Contributed by William A. Goddard III, February 6, 2009 (received for review November 25, 2008)
Abstract
We develop and validate a density functional, XYG3, based on the adiabatic connection formalism and the Görling–Levy couplingconstant perturbation expansion to the second order (PT2). XYG3 is a doubly hybrid functional, containing 3 mixing parameters. It has a nonlocal orbitaldependent component in the exchange term (exact exchange) plus information about the unoccupied Kohn–Sham orbitals in the correlation part (PT2 double excitation). XYG3 is remarkably accurate for thermochemistry, reaction barrier heights, and nonbond interactions of main group molecules. In addition, the accuracy remains nearly constant with system size.
 Becke 3parameter hybrid functional combined with Lee–Yang–Parr correlation functional
 density functional theory
 generalized gradient approximation
 local density approximation
 mean absolute deviation
Density functional theory (DFT) has revolutionized the role of theory by providing accurate firstprinciples predictions of critical properties for applications in physics, chemistry, biology, and materials science (1). Despite dramatic successes, there remain serious deficiencies, for example, in describing weak interactions (London dispersion), which are so important to the packing of molecules into solids, the binding of drug molecules to proteins, and the magnitude of reaction barriers. We propose here a DFT functional that dramatically improves the accuracy for these properties by including the role of the virtual (unoccupied) states.
Solution of the Schrödinger equation leads to the wavefunction, ψ(r_{1}, r_{2}, …, r_{N}) (2), which depends on the 3N space coordinates and N spin coordinates of Nelectrons in the system. Solving for such a wavefunction usually starts with the Hartree–Fock (HF) mean field description involving N selfconsistent 1particle spinorbitals (in a Slater determinant), which is then used as the basis for expanding the wavefunction in a hierarchy of excited Nelectron configurations, by using methods referred to as Møller–Plesset theory (e.g., MP2, MP3, MP4), couplecluster theory (e.g., CCSD(T)), and quadratic configuration interaction theory (e.g., QCISD(T)), etc. These methods are ab initio but suffer from problems of slow convergence with the size of the basis sets and the configuration expansion lengths, preventing scaling to large systems.
In contrast, DFT is formulated in terms of the 1particle density, ρ(r), depending on only 3 spatial coordinates rather than 3N, as the fundamental quantity (3, 4). This dramatically simplifies the process of calculating the structures and properties. However, the exact form of the functional, whose solution will lead to the correct density, is not known. Even so, there has been an evolution of successively better approximations to this functional, that has already provided quite good accuracy for many problems (5–15).
Perdew (16) has formulated the hierarchy of DFT approximations as a “Jacob's ladder” rising from the “earth of Hartree” to the “heaven of chemical accuracy.” The first rung of this ladder is the local (spin) density approximation [LDA, e.g., SVWN (4, 5)] and the second rung is the generalized gradient approximation [GGA, e.g., BLYP (6, 7) and PBE (8)]. Although LDA uses densities ρ(r) as local ingredients, GGA employs both the local densities and their gradients ▿ρ(r). The third rung is termed metaGGA [e.g., TPSS (9)], which expands GGA to include further the kinetic energy density τ, and/or the Laplacian of the density ▿^{2}ρ(r). Up to this third rung, they are all local and multiplicative.
The fourth rung of DFT is a hybrid that introduces nonlocality by replacing some portion of the local exchange energy density with the exact (HFlike) exchange energy density. The most popular such hyperGGA flavor is B3LYP (5–7, 10), which has been shown to provide accurate predictions for thermochemistry of small covalent systems (11). However, B3LYP is poor for the predictions of noncovalent bonding interactions (15) and reaction barrier heights (14), with performance degrading dramatically as system sizes increase (12, 13).
The final fifth rung of Jacob's ladder utilizes the unoccupied Kohn–Sham (KS) orbitals (16) in addition to the occupied KS orbitals. This final rung is expected to allow the heaven of chemical accuracy to be achieved for broad applications. However, no such functional based on first principles (17) and practical for general use has been proposed. Empirical versions (18) have led to promising results for thermochemistry and reaction barriers, but they still fail to account for van der Waals interactions.
Here, we develop a fifthrung functional that incorporates information about the unoccupied KS orbitals [based on the Görling–Levy couplingconstant perturbation expansion to the second order (19)], along with 3 empirical mixing parameters. We demonstrate that this functional is highly accurate for thermochemistry, reaction barriers, and nonbond interactions.
Theory
DFT was placed on a firm theoretical footing by the Hohenberg–Kohn (HK) theorems (3). These HK theorems prove that there exists a total energy functional E[ρ], from which one can obtain the ground state electron density ρ_{0} by minimizing E[ρ] with respect to the density ρ, where ρ_{0} contains all information that can be known about the electronic structure of the system. However, the HK theorems do not specify this true total energy functional.
The most popular implementation of DFT is through the KS method (4), which assumes a noninteracting Nelectron system having the same density as the original manybody system. The KS wavefunction can be expressed exactly as a Slater determinant leading to an exact form for the kinetic energy T_{s} of the noninteracting system and the classic Coulomb energy U. The total energy is then expressed as where V_{ext} is the external potential energy, and E_{xc} is the exchangecorrelation energy, which remains unknown.
The adiabatic connection formalism (10, 20–25) provides a rigorous way to define E_{xc}. It assumes an adiabatic path between the fictitious noninteracting KS system (λ = 0) and the physical system (λ = 1) while holding the electron density ρ fixed at its physical λ = 1 value for all λ of a family of partially interacting Nelectron systems: U_{xc,λ} is the potential energy of exchange correlation at intermediate coupling strength λ. The only problem is that the exact integrand U_{xc,λ} is unknown.
Becke first used this formalism as a practical tool for functional construction (10, 23) by assuming a linear model (23) and taking U_{xc,λ=0} = E_{x}^{exact}, the exact exchange of the KS orbitals, and approximating U_{xc,λ=0} ≈ U_{xc,λ=0}^{LDA}. Becke's halfandhalf functional (23) may be approximated by where we have partitioned E_{xc}^{LDA} = E_{x}^{LDA} + E_{c}^{LDA} and set The popular Becke's 3parameter functional modifies Eq. 5 empirically to obtain Eq. 7 (10): where ΔE_{x}^{GGA} is the gradientcontaining correction terms to the LDA exchange and ΔE_{c}^{GGA} is the gradientcontaining correction to the LDA correlation, whereas {c_{1},c_{2},c_{3}} are constants fitted against selected experimental thermochemical data. The success of Eq. 7 in achieving high accuracy demonstrates that errors of E_{xc}^{DFT} for covalent bonding arise principally from the λ → 0 or exchange limit, making it important to introduce some portion of exact exchange (10, 23–25).
An alternative to fixing the {a,b} parameters in Eq. 4 is to use the Görling–Levy theory of couplingconstant perturbation expansion (19), in which the initial slope (U′_{xc,λ=0}) is defined by the secondorder correlation energy: We may define E_{c}^{GL2} as (19): where ν̂_{ee} is the electron–electron repulsion operator, ν̂_{x} is the local exchange operator, and f̂ is the Focklike, nonlocal exchange operator. We may calculate E_{c}^{GL2} from the KS orbitals with eigenvalues ε, where the subscripts (i, j) and (α, β) denote the occupied and unoccupied KS orbitals, respectively.
Combining Eq. 8 with Eq. 4 leads to: Eqs. 6 and 10 lead to 2 choices of b, which we combine using empirical parameters, {b_{1},b_{2}}, to optimize the functional performance: In principle, E_{c}^{DFT} ≈ (E_{xc}^{DFT} − E_{x}^{exact}) contains a complete description of correlation effects, so that the second term of Eq. 11 may be interpreted as a way to extrapolate the secondorder perturbation to infinite order. Hence, we propose to use an empirical formula of the form: In comparison with the Becke 3parameter scheme (10) of Eq. 7, Eq. 12 is a doubly hybrid DFT that mixes some exact exchange into E_{x}^{DFT} while also introducing a certain portion of E_{c}^{PT2} into E_{c}^{DFT}. Here, E_{c}^{PT2}contains the doubleexcitation contributions of E_{c}^{GL2} (i.e., the first term in Eq. 9). The singleexcitation contributions in E_{c}^{GL2}may not be zero, but we absorb them into E_{c}^{DFT} and in the fitting parameters. Eq. 12 presents a fifthrung functional (R5) that embodies information from both the occupied and the unoccupied KS orbitals as shown in Eq. 9.
In our current applications to test this functional, we calculate the B3LYP wavefunction and use the B3LYP orbitals as the KS orbitals to generate the density and to evaluate the PT2 term. Instead, the original GL2 perturbation theory (19) uses KS orbitals generated from a local exchangecorrelation potential (see Eq. 9). Ref. 26 has shown that B3LYP densities are similar to those from CCSD(T) ab initio wavefunctions (for the molecules discussed in ref. 26). Nevertheless, the eigenvalues from B3LYP, whose potential is nonlocal, might differ considerably from those of the KS orbitals obtained from a local potential. Thus, it could be better to use some different set of KS orbitals.
Here, we adopt the LYP correlation functional but constrain c_{4} = (1 − c_{3}) in Eq. 12. This constraint is not necessary, but it eliminates 1 fitting parameter while excluding compensation from the LDA correlation term. The final 3 parameters {c_{1},c_{2},c_{3}} are determined empirically by fitting only to thermochemical data of the G3/99 set, leading to: We denote this generalized 3parameter functional as XYG3.
Results and Discussion
Heats of Formation (Thermochemistry).
The Gn paradigm developed by Pople and coworkers provides a hierarchy for extrapolating levels of correlation and basis sets to obtain increasingly accurate thermochemistry (11, 12, 27). To adjust the empirical constants in Gn, they developed a database (DB) of accurate experimental heats of formation that are valuable for developing functionals to describe covalent bonding in the main group molecules. In particular, we use the G3 DB of 223 molecules collected in 1999 (the G3/99 set) (12).
Using XYG3 with the 6311+G(3df,2p) basis set to calculate the heats of formation of the G3/99 set leads to a mean absolute deviation (MAD) of 1.81 kcal/mol, substantially better than any other DFT methods (Table 1). For comparison, B3LYP leads to MAD = 4.74 kcal/mol. Indeed the G3 method gives MAD = 1.05, whereas G2 gives MAD = 1.88 kcal/mol but at far higher computational cost.
A recent important development in DFT is the M06 family of functionals (M06, M062X, M06HF, and M06L) (14, 30), which currently provides the highest accuracy with a broad applicability for chemistry. M06, M062X, M06HF are hybrid methods, whereas M06L is a pure DFT. For the G3/99 set, these methods lead to MAD = 4.17 kcal/mol for M06, 2.93 for M062x, and 5.82 for M06L.
B2PLYP is also a doubly hybrid functional that incorporates a perturbation correction as in Eq. 12, but with different parameters of {c_{1} = 0.53, c_{2} = 0.47, c_{3} = 0.27} (18, 31). The first 2 parameters for the exchange part are normalized to 1.0, which reduces the number of independent fitting parameters to 2. The salient difference between B2PLYP and XYG3 is that B2PLYP employs the DFT portion of Eq. 12 to generate the density used to calculate the DFT energy and orbitals from which the PT2 correction is computed. Such a truncated DFT may give density and orbitals that are dramatically different from the real ones. Thus, using just the DFT portion of B2PLYP leads to MAD = 174.2 kcal/mol for the G3/99, whereas the complete B2PLYP method leads to MAD = 4.63 kcal/mol (with our present basis set). Because the doubly hybrid functionals are rooted within the adiabatic connection theorem (20, 21) and the Görling–Levy theory of couplingconstant perturbation expansion (19), we consider it very important to have accurate KS orbitals to provide an accurate density and the zeroorder approximation for perturbation theory.
The G3/99 set consists of 3 subsets of molecules: G2–1 with 56 molecules having up to 3 heavy atoms, G2–2 with 92 additional molecules up to 6 heavy atoms, and G3–3 with 75 additional molecules up to 10 heavy atoms. B3LYP leads to errors that increase dramatically with size (12, 13), with MAD = 2.12 kcal/mol (G21), 3.69 (G22), and 8.97 (G33). B2PLYP [at the 6311+G(3df,2p) level] does not improve over B3LYP, leading to MADs of 1.85 (G21), 3.70 (G22) and 7.83 kcal/mol (G33). M06L gives MADs of 3.76 (G21), 5.71 (G22) and 7.50 kcal/mol (G33). This is significantly improved by M062X, which includes a doubled portion of exact exchange, leading to MADs of 1.89 (G21), 3.22 (G22), and 3.36 (G33) kcal/mol. For XYG3, we obtain MADs of 1.52 (G21), 1.79 (G22), and 2.06 kcal/mol (G33), which exhibits the best description for larger molecules.
Reaction Barrier Heights.
Zhao and Truhlar compiled several benchmark DBs of barrier heights in 2004 (14, 15, 33), including forward and reverse barrier heights for 19 hydrogen transfer (HT) reactions, 6 heavyatom transfer (HAT) reactions, 8 nucleophilic substitution (NS) reactions, and 5 unimolecular and association (UM) reactions. We used the 6311+G(3df,2p) basis set to compute the barriers (see Table 2). Geometries and reference energies were taken from the Truhlar DB web site (14, 15, 33).
DFT methods usually underestimate reaction barrier heights. Table 2 shows MAD errors (kcal/mol) of 14.88 (LDA), 8.71 (PBE), and 4.28 (B3LYP). B2PLYP (MAD = 1.94) leads to a substantial improvement, but M062X (MAD = 1.20) and XYG3 (MAD = 1.02) are remarkably accurate for all types of reactions for a total of 76 barrier heights. This accuracy is comparable with that of the QCISD(T) ab initio method with the same basis set (1.10 kcal/mol). We emphasize that barrier heights are not included in the XYG3 training set but are included in the M06 training set. Probably it is the presence of ≈80% exact exchange in XYG3 that decreases the selfinteraction errors (SIE) of local DFT functionals (25), whereas SIE make local DFT functionals problematic for the stretched partially broken bonds, characteristic of the transition states for chemical reactions.
Accurate potential energy surfaces (PES) are essential for using theory to predict chemical processes, but the accuracy depends critically on the level of the theory. Here, we test various methods for describing the H + CH_{4} → H_{2} + CH_{3} reaction. Because of its important roles in CH_{4}/O_{2} combustion chemistry, this reaction has long been the subject of both experimental and theoretical interest (34). Fig. 1 presents a pointtopoint comparison among the results of various methods along the reaction coordinate. We expect that the CCSD(T) curve should be the most accurate, leading to a barrier of 15.03 kcal/mol. Remarkably, XYG3 predicts the barrier of 15.08 kcal/mol, and is within 0.6 kcal/mol of the CCSD(T) results for the entire reaction path.
The LDA (SVWN) reaction path is qualitatively wrong, predicting a barrierless reverse reaction. HF overestimates the barrier height by 10.89 kcal/mol, whereas HF_PT2, which uses exact exchange plus PT2 correlation, overestimates the endothermicity of the reaction by 6.30 kcal/mol. Here, the tendency that BLYP underestimates the barrier heights is seen clearly in Fig. 1, whereas B3LYP, with inclusion of some exact exchange, leads to improved results, but remains inadequate for PES calculations.
Noncovalent Interactions.
The noncovalent interaction DB from Zhao and Truhlar (14, 15) (NCIE31/05) consists of 6 HB complexes, 7 chargetransfer (CT) complexes, 6 dipole interaction (DI) complexes, 7 weak interaction (WI) complexes, and 5 π–π stacking (PPS) complexes. We tested the performance of the XYG3 functional for these noncovalent interactions using the 6311+G(3df,2p) basis set, with geometries and reference energies taken from the Truhlar DB web site (14, 15).
The errors are summarized in Table 3. We did not include basis set superposition error corrections, which may increase the calculated interaction energies slightly.
We expect that the QCISD(T) ab initio method provides the highest accuracy, leading to a MAD = 0.57 kcal/mol. Remarkably M062X (MAD = 0.30) and XYG3 (MAD = 0.32) have half this error including WI and PPS for which London dispersion dominates. Also M06 (MAD = 0.43) and M06L (MAD = 0.58) perform well for all 5 sets. Note that these nonbond interactions were not included in the XYG3 training set but were included in the M06 training set.
WI and PPS lead to notorious failures for common DFT methods because dispersion interactions are lacking in the correlation functionals. The PT2 term in XYG3 reduces this error, but B2PLYP was not able to describe the PPS complexes. It was suggested that this might be because the PT portion (≈27%) is too small to overcome the repulsion from the DFT part (35). However, we suspect that it is because the orbitals from the truncated DFT in B2PLYP stray too far from the real KS orbitals.
Fig. 2 compares the intermolecular potentials of the CH_{4}–C_{6}H_{6} complexes calculated by XYG3 and CCSD(T) (36), along with some other DFT results. Proper description of such potential energy curves is very important for describing the binding of ligands to biological systems, because steric constraints might prevent the ligand from adopting its optimum geometry. We find that XYG3 results compare extremely well with those of CCSD(T), with deviations generally <0.1 kcal/mol.
Fig. 2B shows the exchange contributions to the noncovalent interaction energy. Here, we see that XYG3 agrees closely with HF, which we expect to be the most accurate. We note here that Slater exchange (S) leads to a spurious well, whereas the GGA correction (Becke88) causes the potential curve to be far too repulsive.
Fig. 2C shows that correlation (attractive) contributions to the noncovalent interaction energy. Here, we see that XYG3 agrees closely with CCSD(T), which we expect to be the most accurate. Note that PT2 by itself is too attractive. Thus, it is the combination of PT2 with LYP that provides the excellent correlation description in XYG3.
That the exchange and correlation parts of XYG3 independently fit what are expected to be the most accurate descriptions indicates that XYG3 gets the right answer for the right reason with a correct description of the fundamental physics.
Summary
We develop here the extension of DFT to the fifth rung of the Perdew Jacob ladder hierarchy. This is done through a construction based on the adiabatic connection formalism using the Görling–Levy couplingconstant perturbation expansion to the second order. This leads to a doubly hybrid functional, XYG3, that uses exact exchange to improve the quality of DFT exchange at the exchange limit (λ = 0), while using both occupied and unoccupied KS orbitals through doubleexcitation contributions from the PT2 term. In this work, we use the KS orbitals and eigenvalues from a selfconsistent B3LYP calculation to compute the PT2 term, which is supplemented with a fraction of LYP to provide the complete correlation energy. Other choices of the virtual orbitals would be possible.
XYG3 contains 3 empirical parameters: (i) the proportion of exact exchange (normalized with the portion of LDA exchange), (ii) the proportion of GGA exchange correction, and (iii) the proportion of PT2 (normalized with the portion of LYP correlation), which determined by using only thermochemical data. In addition to the accuracy of XYG3 for thermochemistry, we find that it is remarkably accurate for the energy along the reaction pathway including reaction barrier heights and for nonbond interactions, neither of which were included in the training set. This suggests that XYG3 captures a consistent description of the physics.
XYG3 does have limitations. Approximate functionals may break the variational principle, leading to energies lower than exact. This can be a serious problem for the PT2 term when there is neardegeneracy of the orbitals as, for example, in the system containing transition metals.
It is also important to consider the scaling of such DFT methods to judge the practicality for application to large systems. Thus, pure DFT methods scale most favorably with size. Including the exact exchange as in B3LYP leads a formal scaling as N^{4}, whereas including the PT2 term leads to a formal scaling as N^{5}, just as for MP2. Linear scaling methods have been developed for MP2 (37–39) that dramatically accelerate calculations for large molecules, and we expect that these can be used with XYG3. See SI for additional information
Acknowledgments
We thank Prof. D. H. Zhang (Dalian Institute of Chemical Physics, Dalian, China) for providing the CCSD(T) results of the potential energy curves for the H + CH_{4} → H_{2} + CH_{3} reaction. This work was supported by National Natural Science Foundation of China Grants 20525311, 20533030, 20423002, and 10774126; Ministry of Science and Technology of China Grants 2007CB815206 and 2004CB719902, with partial support by National Science Foundation Grants ECS0609128 and CTS0608889) and Office of Naval Research (ONR)–Defense Advanced Research Projects Agency Grants PROM N000140610938 and N000140510778). The computation facilities of the Materials and Process Simulation Center (MSC) used in these studies have been supported by grants from the Army Research Office–Defense University Research Instrumentation Program (DURIP) and ONR–DURIP.
Footnotes
 ^{1}To whom correspondence may be addressed. Email: xinxu{at}xmu.edu.cn or wag{at}wag.caltech.edu

Author contributions: X.X. and W.A.G. designed research; Y.Z. and X.X. performed research; Y.Z., X.X., and W.A.G. analyzed data; and Y.Z., X.X., and W.A.G. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0901093106/DCSupplemental.
References
 ↵
 Chong DP
 ↵
 Szabo A,
 Ostlund NS
 ↵
 Hohenberg P,
 Kohn W
 ↵
 Kohn W,
 Sham JL
 ↵
 Vosko SH,
 Wilk L,
 Nusair M
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Perdew JP,
 et al.
 ↵
 Constantin LA,
 Pitarke JM,
 Dobson JF,
 GarciaLekue A,
 Perdew JP
 ↵
 Grimme S
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 MoriSanchez P,
 Cohen AJ,
 Yang WT
 ↵
 ↵
 ↵
 ↵
 ↵
 Zhao Y,
 Truhlar DG
 ↵
 ↵
 ↵
 Zhang LL,
 Lu YP,
 Lee SY,
 Zhang DH
 ↵
 ↵
 ↵
 Kobayashi M,
 Imamura Y,
 Nakai H
 ↵
 ↵
Citation Manager Formats
Sign up for Article Alerts
Jump to section
You May Also be Interested in
More Articles of This Classification
Physical Sciences
Chemistry
Biological Sciences
Related Content
 No related articles found.
Cited by...
 First principlesbased multiscale atomistic methods for input into first principles nonequilibrium transport across interfaces
 Doubly hybrid density functionals that correctly describe both density and energy for atoms
 A fast doubly hybrid density functional method close to chemical accuracy using a local opposite spin ansatz