## 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

# A fast doubly hybrid density functional method close to chemical accuracy using a local opposite spin ansatz

Contributed by William A. Goddard III, September 21, 2011 (sent for review April 15, 2011)

## Abstract

We develop and validate the XYGJ-OS functional, based on the adiabatic connection formalism and Görling-Levy perturbation theory to second order and using the opposite-spin (OS) ansatz combined with locality of electron correlation. XYGJ-OS with local implementation scales as N^{3} with an overall accuracy of 1.28 kcal/mol for thermochemistry, bond dissociation energies, reaction barrier heights, and nonbonded interactions, comparable to that of 1.06 kcal/mol for the accurate coupled-cluster based G3 method (scales as N^{7}) and much better than many popular density functional theory methods: B3LYP (4.98), PBE0 (4.36), and PBE (12.10).

Obtaining chemical accuracy (∼1 kcal/mol) to quantify key chemical quantities (e.g., heats of formation, bond dissociation energies, and reaction barrier heights) using quantum mechanics (QM) has been a major focus in the development of the theory. This has led to, for example, the Gn method (1, 2) that approaches this chemical accuracy. Because G3 is a coupled-cluster based method, it scales on the order of N^{7}, where N measures the system size, limiting to fairly small species for routine use.

The desire to predict unique physiochemical phenomena (e.g., solvation, catalysis, self-assembly, and drug design) in practical (large) systems has brought about a second major focus of theoretical development, leading, for example, to divide-and-conquer formulations to attain more efficient scaling (3), but with much lower accuracy than Gn.

Density functional theory (DFT) in the framework of Kohn-Sham (KS) scheme (4, 5) provides a “shortcut” to the many-body problem. Many density functional approximations (DFAs), provide typical scaling of N^{3} ∼ N^{4}, while yielding significantly more accurate results than Hartree-Fock (HF) theory, the lowest level wave function based method with similar scaling, but they still lead to significant errors for some systems. For example, current DFAs lead to a poor description of London dispersion (van der Waals attraction), which is essential to predict the packing of molecules into solids, and the binding of drug molecules to proteins. These DFAs are also poor in predicting the magnitude of reaction barriers.

In this article, we develop and present a unique functional, XYGJ-OS, that provides a good combination of high accuracy and speed. XYGJ-OS involves a doubly hybrid density functional (DHDF), containing both a nonlocal orbital-dependent component in the exchange term (HF-like exchange), and also information about the unoccupied KS orbitals in the electron correlation part (PT2, perturbation theory up to second order) using the opposite-spin (OS) ansatz to include the locality of electron correlation. XYGJ-OS provides accuracy comparable to that of Gn for the test datasets and speed with N^{2} ∼ N^{3} for the local implementation. Hence XYGJ-OS is both accurate and fast.

## Theory

The Holy Grail in KS-DFT is to find the exact, yet unknown, exchange-correlation functional *E*_{xc}[*ρ*] using density *ρ* as the basic variable (4, 5). In practice, an approximate *E*_{xc} must be adopted, which is often partitioned into the exchange and correlation parts [1] has been extended to include a portion of nonlocal , where the superscript “HF” emphasizes that the exchange part has the same form as in Hartree-Fock theory (6). The exchange-correlation potential *υ*_{xc} in the nth cycle of the self-consistent-field (SCF) process to solve the KS equation is obtained as a functional of *ρ* of the previous cycle.

On the other hand, Görling and Levy (GL) (7, 8) argued that the same KS scheme should work as well in terms of KS orbitals *φ*_{i} and eigenvalues *ε*_{i}. GL proposed a formally exact KS scheme based on perturbation theory, where *E*_{xc} was expressed as: [2][3]Here stands for the GL perturbation theory up to the second order, {Φ_{s,0},*E*_{s,0}} and {Φ_{s,i},*E*_{s,i}} are the wave function and energy for the ground state and the *i*th excited state of an *N*-electron KS system, respectively; is the operator of electron-electron repulsion, and the *υ*_{1} potential may be determined from the “exchange-only” KS equation (8). With knowledge of the potentials *υ*_{j}(*r*), Eq. **2** gives the formally exact exchange-correlation energy as functional of the KS orbitals *φ*_{i} and eigenvalues *ε*_{i}, with which the approximate defined in [**1**] can be compared. However, in practice, such a procedure is difficult to apply to higher than second order due to the unfavorable scaling with the system size, and in many cases perturbation theory is nonconvergent. Thus the scheme must be simplified to make it applicable for including higher-order contributions.

The adiabatic connection method [(ACM), ref. 9, 10], which defines a family of partially interacting N-electron systems with a coupling constant for a fixed *ρ*, provides a powerful tool for developing and understanding *E*_{xc} (6, 11–13). ACM suggests that Eq. **2** up to second order, while being more appropriate for coupling close to zero, is only exact if the potential energy of exchange-correlation depends linearly on the coupling strength, whereas [**1**] is more associated with full electron-electron coupling (13, also see *SI Text* for more discussion). To go beyond the linear approximation and to take advantages of both [**1**] and Eq. **2**, we define an empirical DHDF that combines Eq. **2** (*j* = 2) with [**1**]: [4]Here are the exchange and correlation components within the local density approximation (LDA), and are the corresponding correction terms to LDA within the generalized gradient approximation (GGA). The meta-GGA functionals that include kinetic energy density or the Laplacian of density can also be used in place of GGAs. Such a combination suggests that only the energetically most important double excitation PT2 term in Eq. **2** need be calculated explicitly, the higher-order correlations are embedded into the parameterized terms such as (12). Other examples of DHDFs include MC3BB (14), B2PLYP (15), B2GP-PLYP (16), and ωB97X-2 (17), which are derived and constructed differently (see below and *SI Text* for more discussion).

Several key issues distinguish our approach from GL perturbation theory (7, 8) and other DHDFs (14–23).

Instead of seeking for {

*φ*_{i},*ε*_{i},*υ*_{j}} self-consistently and order by order via Eq.**2**, we use a conventional DFA defined in [**1**] to generate {*φ*_{i},*ε*_{i}} and to evaluate*E*_{xc}defined in Eq.**4**in a post-SCF manner. In fact, the GL perturbation theory, Eq.**2**, was built on KS orbitals with a local effective potential, while we have extended it to include hybrid functionals with generalized KS orbitals (12, 13).differs from in that the singles contribution is not explicitly calculated, but we argue (12, 13) that it can be reasonably absorbed into the parameterization in Eq.

**4**. Higher-order contributions*E*_{c,j}(*j*> 2) are also implicitly included via the parameterization.Neglecting the nonlocal exchange-correlation effects in Eq.

**4**leads to a pure DFA, while neglecting the local exchange-correlation effects, Eq.**4**gives an approximation to GL2. Hence our DHDF may be regarded as an interpolation approach between a pure DFA and GL2, while both of them are in the framework of the KS scheme. As different orbitals are used to construct the PT2 term, neglecting the local exchange-correlation effects in other DHDFs (14–23) will bring back the conventional MP2, such that these DHDFs may be regarded as an interpolation approach between a DFA and MP2, the wave function based lowest level correlation method. These functionals go beyond the framework of the KS scheme (23).

Assuming LDA and GGA in Eq. **4** are SVWN (24, 25) and BLYP (26, 27), respectively, we previously developed XYG3 (12), a DHDF which is remarkably accurate for a wide range of systems and important chemical properties (12, 13, 28–31). Nevertheless, the PT2 term in XYG3 and other DHDFs is evaluated in a way similar to MP2: [5]where the subscripts (*i*, *j*) and (*a*, *b*) denote the occupied and unoccupied KS spin-orbitals, respectively, leading to a formal scaling as N^{5}, as opposed to a formal scaling of N^{4} as in B3LYP (24–27, 32, 33). This unfavorable scaling raises an issue for the practicality to apply DHDFs to large systems.

Density fitting approximations have often been used in electronic structure theories to reduce computational expense. In the so-called “resolution-of-the-identity” RI-MP2 method (34, 35), the product of occupied and virtual orbitals (*ia* pair) is expanded with auxiliary functions, such that numerous 4-center 2-electron integrals based on molecular orbitals (MO) are replaced by fewer 3 and 2-center integrals with cheaper transformation from atomic orbitals (AO) to MO. RI-MP2 is about 5–20 times faster than conventional MP2 (36, 37). Nevertheless, RI-MP2 reduces only the prefactor, it does not change the scaling. Similarly, we have RI-XYG3, which retains the original accuracy but is faster for small systems with large basis sets.

Here we propose a unique OS ansatz for DHDF, that yields a balanced description of nonlocal correlation effects while considerably reducing computational time. Our OS ansatz is motivated by the observation that the most important electron correlation effects involve correlations of the OS electrons in the same orbital. The OS ansatz leads to N^{4} scaling (36, 37) [using auxiliary basis expansions (34, 35) and Laplace quadrature approximations (38)]. It should be emphasized that as the same-spin (SS) correlation is very important in accurate description of open-shell systems and magnetic properties, such contributions cannot be simply neglected. In XYGJ-OS, the SS correlation effects are included within the standard DFA.

In recognition of the “nearsightedness of electron correlation” as emphasized by Kohn (39), we then build upon XYGJ-OS to introduce the local approximation for the OS electron correlation by utilizing the sparsity of the RI expansion coefficients, integral matrices, and Laplace transform matrices. The local implementation of XYGJ-OS scales as N^{3} while retaining the accuracy of the original XYGJ-OS (see *SI Text* for more discussion).

Thus our proposed functional form (XYGJ-OS) is [6]In Eq. **6** we normalize the HF exchange and Slater exchange (24), while eliminating the contribution. The correlation part consists of (25), (27), and , where the first term includes both the SS and OS effects while the second and third terms include only OS components. Our concept is that the combination of VWN, LYP, and PT2-OS yields a balanced description of both local and nonlocal spin dependent correlation terms. To determine the optimal four parameters in Eq. **6**, we use the experimental heats of formation (HOF) data for the G3/99 set of 223 molecules (1, 2) as the training set, leading to {*e*_{x},*e*_{VWN},*e*_{LYP},*e*_{PT2}} = {0.7731,0.2309,0.2754,0.4364}.

## Results and Discussion

We emphasize that we use the fully optimized B3LYP orbitals to generate the density and to calculate each term in Eqs. **4** or **6** (12, 13, 29–31). But the choice of LDA = SVWN and GGA = BLYP in Eqs. **4** or **6**, as well as using B3LYP orbitals as input, is not unique. We find that any conventional DFAs defined in [**1**] can also serve the same purpose, leading to similar performance, albeit with a reoptimized set of mixing parameters. In our method, we assume that our KS wave function is the zeroth-order wave function in the GL perturbation theory that gives the correct ground-state electron density. Instead, MC3BB is a multicoefficient method, which mixes the total energies from the conventional MP2 and a conventional DFA calculation. Hence there are two independent SCF calculations in the MC3BB type of DHDFs, which lead to two sets of different orbitals, yielding two different densities. In MC3BB, the SCF-HF orbitals are used for MP2 evaluation (14), while it is well known that HF wave function is the one-determinantal wave function which gives the lowest expectation value with the fully interacting Hamiltonian. It has been shown recently (23) that the B2PLYP family of functionals works in a way which is technically similar to MP2. These functionals adopt orbitals that minimize the one-determinantal wave function based on the so-called density-scaled one parameter hybrid approximation. Hence, the density from the B2PLYP type of functionals is by construction not meant to be the true density. In contrast to what was originally proposed (15), there is no singles’ contribution by construction in the B2PLYP family of functionals, and its theoretical basis is provided by the multideterminant extension of the Kohn-Sham scheme (23). On the other hand, the key idea of the XYG3 type of functionals is to combine the GL perturbation theory and the standard KS scheme in the framework of ACM. In XYG3, only the energetically most important double excitation PT2 term in Eq. **2** is calculated explicitly using orbitals generated from a conventional DFA in [**1**]. Our present OS ansatz further reduces the computational cost by only calculating the most important electron correlation effects contributed by the OS electrons in the same orbital. XYGJ-OS presents a unique combination of speed and accuracy. Hence, we propose that the current DHDFs shall be categorized into three types (13), as represented by MC3BB (14), B2PLYP (15), and XYG3 (12). The former two go beyond the KS scheme.

It is possible to carry out a SCF calculation for any orbital-dependent *E*_{xc} using, for example, the direct optimization approach to compute the optimized effective potential (OEP) as proposed by Yang and coworkers (40). It was found that such a SCF procedure can often lead to unphysically unbound state with too low total energy mainly due to the near-degeneracy of the orbitals in the single excitation terms. This unbound issue can be largely remedied by calculating only the double excitation terms in the post-SCF manner (40). This procedure is indeed the calculation approach adopted by the XYG3 type of functionals. Furthermore, we note that extra terms will appear when formulating the analytical gradients for doing geometry optimization with our functionals, which, however, will not impose any practical difficulty in implementation as compared to other types of DHDFs.

We find that XYGJ-OS is remarkably accurate for a broad range of systems and important chemical properties (1, 2,41–45) other than HOF which are *not used in fitting the parameters*. Table 1 (more details are in *SI Text*) compares the overall performance of some representative DFT methods, showing that XYGJ-OS is the best or nearly best for essentially all properties, leading to chemical accuracy (1.28 kcal/mol) comparable to the G3 theory (which contains four empirical parameters involving the number of electron pairs) (1.06) and much better than MP2 (7.49 kcal/mol).

### Heats of Formation.

The 223 molecule G3/99 set (1, 2) provides a test of accuracy for the main group covalent systems. XYGJ-OS gives mean absolute deviation (MAD) of 1.65 kcal/mol, lying between those of G2 (1.89) and G3 (1.06) theories. Note that the MADs for HOF listed in Table 1 associated with XYG3, B2PLYP, B2PLYP-D, and ωB97X-2(LP) are taken from the original refs. 12, 17, 22, which were obtained by using the way these functionals were parameterized. It was shown that HOF calculations with DHDFs are most prone to the basis set effects (29). The results with the G3Large basis set (1) used for optimization of the XYGJ-OS functional can be found in *SI Text*, where severe basis set dependences are clearly shown.

### Charged Species.

Charged species are not included in our training set. The G2-1 set (2) for ionization potential (IP) contains 14 atoms and 24 molecules. XYGJ-OS gives a MAD of 1.23 kcal/mol, being one of the best DFT methods for calculating IP. Over the 25 cases in the G2-1 set for electron affinity (EA) (2), XYGJ-OS leads to MAD = 1.97 kcal/mol. Generally, increasing the size of the basis set will increase the accuracy for EA calculation with DHDFs. For the eight systems for proton affinity in the G2 and G3 sets (1), XYGJ-OS leads to MAD = 1.68 kcal/mol, comparable to the performance of conventional DFAs.

### Bond Dissociation Enthalpy (BDE).

BDE is arguably the most important property for chemistry, carrying additional information other than HOFs (13). XYGJ-OS leads to MAD = 0.71 kcal/mol for the 92 reactions in the BDE92/07 set (40)! This accuracy surpasses the MAD of the most accurate ab initio methods, G2 (1.80) and G3 (1.09 kcal/mol). Indeed XYGJ-OS leads to half the MAD of XYG3 (1.57 kcal/mol), providing strong support for the OS ansatz and the local approximation.

### Reaction Barrier Height (RBH).

We use the NHTBH and HTBH datasets of Zhao and Truhlar as testing sets (42, 43) of RBHs. NHTBH contains six heavy-atom transfer reactions, eight nucleophilic substitution reactions, five unimolecular and association reactions, while HTBH comprises solely 19 hydrogen transfer reactions. XYGJ-OS leads to MAD = 1.03 kcal/mol for the two sets together, a significant improvement over B3LYP (4.55 kcal/mol) and MP2 (4.67 kcal/mol). Note that the dispersion corrected methods [e.g., B3LYP-D3 (46)] may deteriorate the performance for RBH calculations.

We test various methods in describing the whole H + CH_{4} → H_{2} + CH_{3} reaction path using the coupled-cluster method with single and double as well as perturbative triple excitations (CCSD)(T)/6-311++G(3df,2pd) data (43) as the reference. The results are depicted in Fig. 1*A*. XYGJ-OS results are nearly identical to the XYG3 and CCSD(T) results before the barrier. But XYGJ-OS overestimates the reaction endothermicity by 1.21 kcal/mol.

### Nonbonded Interaction (NBI).

We use the NCIE dataset of Zhao and Truhlar (42, 43) to test the XYGJ-OS performance for the description of NBIs. The NCIE set includes six hydrogen bond complexes, seven charge-transfer complexes, six dipole interaction complexes, seven weak interaction complexes, and five *π*-*π* stacking complexes. XYGJ-OS (MAD = 0.35 kcal/mol) does quite well for NBI, including the London dispersion dominant systems in the NCIE dataset. For the S22 set (47), which contains larger molecules (e.g., uracil dimer, adenine-thymine complexes) that are more biologically related, XYGJ-OS (MAD = 0.36 kcal/mol) leads to similar accuracy as the dispersion corrected methods (22, 46). On average, B3LYP-D3 improves B3LYP by 0.34 kcal/mol for the NCIE set. While significant improvements occur for the weak interaction complexes, and the *π*-*π* stacking complexes, B3LYP-D3 significantly overbinds the hydrogen bond complexes, and the charge-transfer complexes (see *SI Text* for more discussion).

We test various methods in describing the intermolecular potentials of the CH_{4}-C_{6}H_{6} complex calculated by various methods. The CCSD(T) results at the complete basis set (CBS) limit are used as reference (44). As shown in Fig. 1*B*, XYGJ-OS data are nearly on top of those of XYG3.

### Scaling.

Equally important is that XYG3 and other DHDFs (15–19) scale formally with the fifth power of the size, while XYGJ-OS scales formally with the fourth power, and for the local implementation it scales formally with the third power. In contrast, the most accurate ab initio methods, CCSD(T), and G3, scale as the seventh power. Indeed, Table 1 shows that the total computational time for *C*_{100} chains and diamondoids with local XYGJ-OS is just 3 to 4 times that for B3LYP (a part of XYGJ-OS), with a scaling of N^{2.1}for n-alkanes (See *SI Text* for more details).

## Summary

We have developed and validated here a unique doubly hybrid density functional, XYGJ-OS, using Walter Kohn’s insight about “nearsightedness” of electron correlation by including explicitly only the correlation between electrons of OS and then only the parts that are in the same region of space. We show that XYGJ-OS achieves nearly chemical accuracy (1.28 kcal/mol) with computational costs scaling as N^{3}. This unique combination of high accuracy and speed leads to a practical level of calculation while attaining chemical accuracy for large molecular systems.

## Acknowledgments

X.X. acknowledges the support of NSF of China (91027044, 21133004, 20923004), and the Ministry of Science and Technology of China (2007CB815206, 2011CB808505). Y.J. acknowledges the support of National Research Foundation (NRF) of Korea (2010-0023018, 2010-0029034, and 2010-0029728) by the Ministry of Education, Science and Technology (MEST) of Korea. W.A.G. acknowledges the support of USA National Science Foundation (NSF) (ECS-0609128, CTS-0608889), the Center for Catalytic Hydrocarbon Functionalization (Department of Energy, Basic Energy Sciences Award DE-SC0001298), and Office of Naval Research-Defense Advanced Research Projects Agency [ONR-DARPA_ (PROM N00014-06-1-0938 and N00014-05-1-0778)]. Y.J. and W.A.G. have also been supported by the World Class University (WCU) (NRF R-31-2008-000-10055-0) program funded by the Ministry of Education, Science and Technology.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. E-mail: xxchem{at}fudan.edu.cn, or ysjn{at}kaist.ac.kr, or wag{at}wag.caltech.edu.

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

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1115123108/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- ↵
- Hohenberg P,
- Kohn W

- ↵
- Kohn W,
- Sham LJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Mori-Sanchez P,
- Cohen AJ,
- Yang WT

- ↵
- Zhang Y,
- Xu X,
- Goddard W A III.

- ↵
- Zhang IY,
- Xu X

- ↵
- ↵
- Grimme S

- ↵
- Karton A,
- Tarnopolsky A,
- Lamère J-F,
- Schatz GC,
- Martin JML

- ↵
- Chai JD,
- Head-Gordon M

- ↵
- Sancho-García JC,
- Pérez-Jiménez AJ

- ↵
- ↵
- ↵
- ↵
- ↵
- Sharkas K,
- Toulouse J,
- Savin A

- ↵
- ↵
- Vosko SH,
- Wilk L,
- Nusair M

- ↵
- ↵
- ↵
- Zhang IY,
- Luo Y,
- Xu X

- ↵
- Zhang IY,
- Luo Y,
- Xu X

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Mori-Sánchez P,
- Wu Q,
- Yang WT

- ↵
- Wu JM,
- Xu X

- ↵
- ↵
- ↵
- Zhang LL,
- Lu YP,
- Lee S-Y,
- Zhang DH

_{4}reaction. J Chem Phys 127:1–7. - ↵
- ↵
- Grimme S,
- Antony J,
- Ehrlich S,
- Krieg H

- ↵

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.