## 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 functionals that correctly describe both density and energy for atoms

Edited by John P. Perdew, Temple University, Philadelphia, PA, and approved January 19, 2018 (received for review July 22, 2017)

## Significance

It was commented that some modern density functional theory functionals may be giving the correct energies for the wrong reason and functionals that produce “right” energies from “wrong” densities fail twice. However, the route to the “heaven” of the “Jacob’s ladder” is incomplete in Medvedev et al. [Medvedev MG, et al. (2017) *Science* 355:49–52], as the authors did not provide any information on the top fifth-rung functionals. We fill the gap to show here that the XYG3 type of doubly hybrid density functionals (xDHs) can correctly describe both density and energy for the same atomic set used by Medvedev et al., such that electronic structure calculations for molecular sciences of the main group elements can be carried out reliably by using xDHs.

## Abstract

Recently, it was argued [Medvedev MG, et al. (2017) *Science* 355:49–52] that the development of density functional approximations (DFAs) is “straying from the path toward the exact functional.” The exact functional should yield both exact energy and density for a system of interest; nevertheless, they found that many heavily fitted functionals for molecular energies actually lead to poor electron densities of atoms. They also observed a trend that, for the nonempirical and few-parameter functionals, densities can be improved as one climbs up the first four rungs of the Jacob’s ladder of DFAs. The XYG3 type of doubly hybrid functionals (xDHs) represents a less-empirical and fewer-parameter functional on the top fifth rung, in which both the Hartree–Fock-like exchange and the second-order perturbative (MP2-like) correlation are hybridized with the low rung functionals. Here, we show that xDHs can well describe both density and energy for the same atomic set of Medvedev et al., showing that the latter trend can well be extended to the top fifth rung.

Density functional theory (DFT; refs. 1⇓–3) is now the leading electronic structure method for atoms, molecules, and extended systems. The bedrock on which DFT stands is the Hohenberg–Kohn theorems (1), which proved that all information for the system of interest is contained in its ground state electron density

Over decades, various approximations to the exchange-correlation energy have been developed (e.g., refs. 3 and 6⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–24). They are constructed according to different philosophies, such that an exchange-correlation functional can be either nonempirical or empirical. The former is formulated only by satisfying some physical and mathematic conditions (3, 6, 7, 10, 11, 13, 16, 17, 24), while the latter is optimized by fitting to the known results of atomic or molecular properties, mostly energy differences (8, 9, 12, 14, 18⇓⇓⇓⇓–23). It is generally believed that the nonempirical approach offers a deeper understanding of the exact functional, whereas the empirical approach tends to do better for systems in the range of fitting than for systems outside this range.

The DFA development toward the “heaven” of chemical accuracy can well be described by the “Jacob’s ladder” of Perdew with five rungs (25). An earlier milestone is the progress from local density approximation [LDA; e.g., Slater (6) and VWN5 (7)] on the first rung to generalized gradient approximation [GGA; e.g., B88 (8) and LYP (9)] on the second rung in the 1980s. The later major advance came in the 1990s with the construction of the nonempirical GGA functional PBE (10) and the inclusion of a fraction of Hartree–Fock-like (HF) exchange, leading to hybrid functionals (14⇓⇓⇓⇓–19) such as B3LYP (14, 15) and PBE0 (16, 17) on the fourth rung. More recent achievements after the 2000s include meta-GGA functionals on the third rung [e.g., TPSS (11), M06-L (12), and SCAN (13)], as well as their hybrid forms [e.g., M06-2X (18)] on the fourth rung. Significant improvements have been well documented for the first four rungs of functionals in predicting accurate energy properties such as heats of formation, bond dissociation enthalpies, reaction barrier heights and complexation energies of nonbonded complexes, etc. (18, 19).

However, Medvedev et al. (26, 27) showed in their recent reports that improvements of the first four rung DFAs in energies are not always on par with their improvements in densities. For a set of isolated atoms and atomic cations with 2, 4, or 10 electrons—i.e., Be^{0}, B^{3+}, B^{+}, C^{4+}, C^{2+}, N^{5+}, N^{3+}, O^{6+}, O^{4+}, F^{7+}, F^{5+}, Ne^{8+}, Ne^{6+}, and Ne^{0}—their results based on 128 DFAs demonstrated that some of the modern functionals on the fourth rung, in particular those based on advanced fitting schemes, produce electron densities that are even less accurate than those from LDAs on the first rung. As energy and density are two sides of one coin, the exact functional shall yield the exact energy of the system from its exact density. Such a mismatch between energy and density suggests a departure from the exact functional (26, 27). This argument adds to the long-standing debate on how DFAs can be improved. Thus, it was commented that some modern empirical functionals may be giving the correct energies for the wrong reason (28, 29), and functionals that produce “right” energies from “wrong” densities fail twice (30).

It has to be noted that the missing information on the fifth rung is very important, as it is the top rung of the Jacob’s ladder, which, in principle, is most close to the heaven. It has been shown that the fifth-rung functionals—e.g., the XYG3 type of doubly hybrid (xDH for short) functionals (21⇓–23)—by properly considering the information of unoccupied orbitals, can produce energies that are often approaching chemical accuracy for systems of different bonding natures for the main group elements (31⇓–33). Herein, we show that densities can also be well produced with the xDH functionals, which are concurrently more accurate than those produced by lower-rung functionals, demonstrating that electronic structure calculations for molecular sciences of the main group elements can be carried out reliably by using xDHs.

## Results

Medvedev et al. (26) have examined the electron density distributions of a set of atomic species, as results for atoms are unambiguous and are sacred to molecular science as uniform electron gases are to extended systems. Hence, a calculation method that fails to describe atoms correctly can have internal problems, which can be expected to propagate to larger systems. To compare the density distributions between different methods, Medvedev et al. (26) have calculated the local electron density (RHO; ρ), the gradient norm of RHO (GRD; |∇ρ|), and the Laplacian of RHO (LR; ∇^{2}ρ). These correspond to the basic ingredients in constructing the LDA, GGA, and meta-GGA functionals, respectively. They adopted the root mean square difference (rmsd) to measure the deviations of the calculated DFT values from the all-electron-coupled cluster singles and doubles, CCSD(full), references between the corresponding radial distribution functions. The quintuple-ς spherical basis sets, aug-cc-pωCV5Z, were used in all their calculations (26, 34). Some of their results are summarized in Table 1, while more can be found in their supporting materials.

Medvedev et al. (26, 27) concluded that (*i*) functionals constructed with little or no empiricism tend to produce more accurate electron densities than highly empirical ones; (*ii*) at the level of little or no empiricism, the accuracy of the density tends to increase along with the complexity of the DFA as one ascends the Jacob’s ladder from LDA to GGA to meta-GGA to hybrids; and (*iii*) recent DFAs on average are shifting away from the exact functional.

To ensure a strict comparison, we follow Medvedev et al. (26), using the same basis set to correlate all electrons and report the normalized rmsd errors on density distributions using the same median rmsd errors for each descriptor (26). For the 128 DFAs considered by Medvedev et al. (26), they found that none of the DFT methods can outperform the MP2 method. The trends from our results for the xDHs (21⇓–23) are immediately clear in Table 1. The xDHs are clearly superior to MP2. In fact, the results from XYGJ-OS (22) and xDH-PBE0 (23) are comparable to those from MP3.

While more details may be found in *SI Materials and Methods*, Table 1 gives the mean normalized rmsd errors for RHO, GRD, and LR, respectively, as well as the corresponding maximum normalized errors for the xDH functionals (21⇓–23). The former shows how well the given method behaves on average as described by the three descriptors of RHO, GRD, and LR; the latter indicates how bad the worst-case behavior is. The net maximum normalized error over all atoms and descriptors on density, which was the indicator used by Medvedev et al. (26) to rank the functional, is marked as MaxD in bold in Table 1. For the mean errors, MP2 gives 0.590, 0.406, and 0.255 for RHO, GRD, and LR, respectively, whereas the corresponding xDH-PBE0 values are 0.388, 0.318, and 0.253, being overall superior. For the maximum errors, MP2 gives 1.414, 1.512, and 1.116 for RHO, GRD, and LR, respectively, whereas the corresponding xDH-PBE0 values are 0.884, 0.928, and 0.747, again being overall superior. The worst case all occurs at Ne^{6+} or Ne, as was also observed for the lower-rung functionals (26). The calculated results are also presented for the averaged normalized errors over all descriptors and atoms for a given number of electrons (n_{e} = 2, 4, or 10). For 2e systems, none of the DFT methods are competitive, whereas even HF is more advantageous. With the increase of electron numbers, DFAs start to win out (35).

Table 1 also gives results from B2PLYP (20), a pioneer and well-established DH functional. As shown by the data in Table 1, B2PLYP on average is nearly as good as MP2 for the same atomic set considered (26). A recent analysis on the correlation of bond energies and bond densities for diatomic molecules also suggested the good behavior of the B2PLYP functional (36).

To relate energy behaviors with density errors directly within the same set that Medvedev et al. (26) have used, Kepp (35) has recently calculated the energy differences for the ion pairs of B^{3+}/B^{+}, C^{4+}/C^{2+}, N^{5+}/N^{3+}, O^{6+}/O^{4+}, F^{7+}/F^{5+}, and Ne^{8+}/Ne^{6+}, all of which correspond to ionization processes from 1s^{2}2s^{2} to 1s^{2}. Taking the B^{3+}/B^{+} pair as an example, the exact energy difference can be obtained from the sum of the second and third ionization potentials (IPs) of boron (35, 37):^{3+}/B^{+} to Ne^{8+}/Ne^{6+}, the experimental IPs can be found in the National Institute of Standards and Technology database (37), which shall provide the reference data. Nonetheless, as the ion charges increase and the electron densities become increasingly compact, the experimental IPs contain increasingly important relativistic effects. Applying the Douglas–Kroll approximation (38) at the HF level, the relativistic correction is estimated to be 0.017 for B^{3+}/B^{+}, which increases to 0.597 eV for Ne^{8+}/Ne^{6+}. Hence, the reference data for the ion-pair energy differences we used are from the experimental IPs corrected by the corresponding calculated relativistic corrections. Such a scalar relativistic correction is particularly effective in the present cases with zero spin and angular momentum (35). This is confirmed by the data in Table 2, which show that nonrelativistic CCSD(full) calculations lead to an rmsd of 0.017 eV with the maximum error of 0.022 eV, in good agreement with the reference data after the scalar relativistic corrections.

Table 2 also summarizes the calculated ion pair energy difference results from some other lower-level ab initio methods. From the data in Table 2, it is clear that there is a systematic improvement from the lowest level, HF, to the highest level, CCSD, used in the present work, either in terms of an individual ion pair, the rmsd from the whole set, or the maximum error. Hence, rmsd decreases from 2.794 (HF), to 1.002 (MP2), to 0.478 (MP3), to 0.219 (MP4SDQ), to 0.017 eV (CCSD).

Errors are actually not at all systematic for the DFT methods. For example, LDA (SVWN5) gives a rmsd of 1.882 eV, which is smaller than those from meta-GGAs (e.g., M06-L) and hybrid functionals [e.g., ωB97X-D (19)]. The latter two give rmsds of 2.649 and 2.912 eV, respectively. Here, we have chosen BLYP (rmsd = 0.329 eV) and PBE (rmsd = 0.569 eV) as the GGA representatives, both of which behave quite satisfactorily for this ion-pair energy difference set. Kepp (35) has examined GGA using BP86 (8, 39) as a representative. The rmsd associated with BP86 was found to be 0.861 eV (35), which is more than doubled compared with the value from BLYP. Table 2 also shows that nonempirical functionals do not necessary perform better than their empirical counterparts. This is clear when one makes the comparison for results between PBE and BLYP (0.569 vs. 0.329 eV) or PBE0 and B3LYP (0.785 vs. 0.307 eV), whereas it is certainly true that TPSS performs better than M06-L (0.726 vs. 2.649 eV). It is also noted that rmsd increases from PBE (0.569) to TPSS (0.726) to SCAN (1.955 eV) for this ion-pair energy difference set, although more constraints have been fulfilled in functional construction along this series.

DH functionals in general behave well. The resultant rmsds are 0.522 eV for B2PLYP and 0.562 and 0.405 eV for XYGJ-OS and XYG3, respectively. The rmsd associated with xDH-PBE0 is only 0.189 eV, with its maximum error of 0.333 eV, which clearly surpasses MP2 and is comparable to MP3 and MP4SDQ.

## Discussion

Fig. 1 summarizes the results which are sorted according to the relative maximum normalized errors on density distributions [i.e., rMaxD = MaxD(Calc.) − MaxD(MP2)] as calculated from the data in Table 1. The MP2 method is used as the reference, which is marked as a star. Therefore, methods degrade from left to right, while a negative or positive rMaxD indicates a method that is better or worse, respectively, than MP2 in terms of density for the Medvedev’s atomic set. When does a functional correctly describe both density and energy? The relative maximum error on the ion-pair energy differences is labeled as rMaxE, in Fig. 1 for each method, which is calculated as rMaxE = MaxE(Calc.) − MaxE(MP2) using data listed in Table 2. The MP2 method is again used as the reference. Hence, a negative or positive rMaxE indicates a method that is better or worse, respectively, than MP2 in terms of energy for the ion-pair set. It is only when rMaxD and rMaxE are both negative that indicates a method better than MP2 for both density and energy. From Fig. 1, we see that all three xDHs fulfill quality requirements, while xDH-PBE0 is the best DFT method examined so far for this atomic set with a concurrent improvement in both density and energy. Note that a positive rMaxD but a negative rMaxE indicates a method in which a poor density yields a good energy.

It is now well recognized that a DFA contains two sources of errors: One is energy-driven, and the other is density-driven. Usually, energy-driven error is dominant (4), which is of the primary significance. In some specific cases, density-driven error becomes significant (40). In fact, these two errors are intercalated in the standard KS scheme: A poor exchange-correlation functional leads to the poor energy, while its functional derivative yields to a poor potential which yields the poor density. Hence, energy and density are related by virtue of the KS potential:*SI Materials and Methods* for more details). Hence, the xDHs ascend and descend the ladder at will for energy and density, having a better chance that both energy and density can be well described.

Previously, we have also shown that the xDH functionals lead to the nearly exact linearity behavior (47) for fractional charges, surpassing other conventional functionals (5, 31, 41, 42). This accounts for the generally improved performances of the xDH functionals in the description of thermochemistry and thermochemical kinetics for the main group chemistry (21⇓⇓–24, 31⇓–33). On the other hand, hybrid functionals on the fourth rung fall short, where the nonlocal HF-like exchange is not balanced by the nonlocal MP2-like correlation as it is for DH functionals on the fifth rung.

The Jacob’s ladder not only helps to categorize the existing approximate functionals, but also is very useful to direct the future for functional development. In this perspective, we would like to conclude that a continuing improvement to the top-rung functionals is the way toward the exact functional. The current construction of DH functionals suffer from the same problem as MP2 for systems with small energy gaps. The other ways of utilizing the information of unoccupied orbitals, e.g., renormalization of the MP2 term (48, 49), random phase approximation instead of MP2 (50), shall help to extend the applicability of the DH functionals.

## Materials and Methods

Doubly hybrid (DH) functionals (20⇓⇓⇓–24) provide one feasible way to utilize the information of unoccupied orbitals, which, in addition to the hybridization of the HF-like exchange

A unique feature for xDHs is that all energy terms associated with the DH functionals are evaluated by using the SCF orbitals from some lower-rung functionals. Specifically, orbitals and densities from B3LYP (14, 15) have been chosen for XYG3 and XYGJ-OS, while those from PBE0 (16, 17) have been chosen for xDH-PBE0. Hence, in the xDH functionals, a lower-rung DFA on the Jacob’s ladder is adopted as an SCF functional for its efficiency and reliability to generate orbitals and densities, while a top-rung DH functional is adopted for final energy evaluation for its improved accuracy due to the introduction of the nonlocal correlation effects. Both the SCF functional and the final energy functional in the xDH functionals can be considered as approximations to the exact, yet-unknown KS functional (31⇓–33). Furthermore, the opposite-spin (OS) ansatz (51) was adopted to simplify the MP2-like term in XYGJ-OS and xDH-PBE0, leading to a favorable O(N^{3} ∼ N^{4}) scaling for large systems.

Even though there are several empirical steps, hybrid functionals such as B3LYP and PBE0 have their roots in the adiabatic connection [AC (52, 53)] formalism. In the same spirit, the xDH functionals have been rationalized and constructed by using the AC formalism and by invoking the Görling–Levy coupling-constant perturbation expansion theory [GL2 (54)]. In terms of the AC-perturbation formalism,

We note that both B3LYP and PBE0 are on their list of the best methods, in the *Science* paper of Medvedev et al. (26), yielding the best densities, while the xDH functionals use these orbitals and densities as the input, which is advantageous by construction. In fact, the B3LYP and PBE0 densities can be updated for an even better accuracy by using a modified Z-vector method (45, 46). The Z vector is obtained from the xDH energy gradient calculations to take into account both the perturbation effects from the DFA part and the MP2-like part, coupling the SCF functional and the DH energy functional together. The Z vector acts as a relaxation correction to the one-particle density matrix obtained from the SCF functional (see *SI Materials and Methods* for more details). The sum of two yields the relaxed xDH density matrices, which are then converted to natural orbitals to compute the corresponding radial distribution functions for ρ, |∇ρ|, and ∇^{2}ρ (55, 56).

B2PLYP (20) is a pioneer DH functional, where a hybrid DFA with partial correlation (e.g., **7**) is used as the SCF functional to generate orbitals and densities, whose total energy is then augmented with a scaled MP2-like term,

The B2PLYP type of DH functionals have been rationalized in terms of the multideterminant extension of the KS scheme with some additional approximations, e.g., on density scaling (57). Note that B2PLYP (20) represents a different type, compared with the xDH functionals, where the **7** is not obliged to give the exact density (41, 42, 57), contrary to the role played by B3LYP or PBE0 in the xDHs. The relaxed densities for B2PLYP can also be calculated via the Z-vector method.

All our calculations were carried out with a local version of the NorthWest Computational Chemistry (NWChem) software package (58).

### Note.

(*i*) When submitting the revised manuscript, we became aware of ref. 59, which found that the DH functionals are slightly less accurate than the MP2 method. Ref. 59 only tested the B2PLYP type of functionals, for which their results agree with the results shown in our present work. We have tested here also the XYG3 type of functionals. The xDH functionals start with a standard hybrid functional (such as B3LYP and PBE0) for SCF whose density is already quite good, and then the xDH energy functionals further update the density to an even better one. Therefore, the final density of an xDH can be better than MP2 as shown here. (*ii*) A very recent work (60) showed that a simpler aug-cc-pV5Z basis set (34) gives very different ranking of various methods. We follow their work and check the xDH performances with aug-cc-pV5Z. We find that the xDH functionals are consistently better than MP2 using both basis sets (*SI Materials and Methods*). (*iii*) There are recent articles discussing some limitations of finite basis set DFT calculations on densities (61, 62). Some preliminary study on the basis set dependences can be found in *SI Materials and Methods*. While the accuracy of the electron density of DFAs is still relatively unstudied, there exist many benchmarking studies (e.g., refs. 18, 19, and 31⇓–33) on the accuracy for molecular energy differences. A recent exhaustive study (63) confirms the remarkable accuracy of the DH functionals, as well as the good performance of the nonempirical lower-rung functionals like SCAN and TPSS in general. (*iv*) Computed densities with xDH functionals for this article may be accessed at http://www.xdft.org/Articles/Supplements/doi_10.1073_pnas.1713047115/.

## Acknowledgments

We thank Dr. Yi Luo, Dr. Dong H. Zhang, and Dr. Xueming Yang for insightful discussion; Dr. Michael G. Medvedev for providing reference data by CCSD and MP2, as well as the updated median errors; Dr. Jianwei Sun for providing the SCAN data; and Dr. Niri Govind and Dr. Álvaro Vázquez Mayagoitia for help with the use of NWChem. This work was supported by National Natural Science Foundation of China Grants 21688102 and 91427301; and The Science Challenge Project Grant TZ2018004.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: xxchem{at}fudan.edu.cn.

Author contributions: X.X. designed research; N.Q.S., Z.Z., and X.X. performed research; X.X. analyzed data; N.Q.S. modified the NWChem package for calculating the relaxed densities of the xDH functionals and generating the corresponding natural orbitals; Z.Z. modified the NWChem package to generate the wfn-formated files for density analyses; and X.X. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

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

Published under the PNAS license.

## References

- ↵
- ↵
- Levy M

- ↵
- ↵
- ↵
- Su NQ,
- Xu X

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Zhao Y,
- Truhlar DG

- ↵
- ↵
- ↵
- Zhang Y,
- Xu X,
- Goddard WA 3rd

- ↵
- Zhang IY,
- Xu X,
- Jung Y,
- Goddard WA 3rd

- ↵
- ↵
- Su NQ,
- Xu X

- ↵
- Perdew JP

- ↵
- Medvedev MG,
- Bushmarinov IS,
- Sun J,
- Perdew JP,
- Lyssenko KA

- ↵
- Medvedev MG,
- Bushmarinov IS,
- Sun J,
- Perdew JP,
- Lyssenko KA

- ↵
- Hammes-Schiffer S

- ↵
- Korth M

- ↵
- Graziano G

- ↵
- Su NQ,
- Xu X

- ↵
- Zhang IY,
- Xu X

- ↵
- Su NQ,
- Xu X

- ↵
- ↵
- Kepp KP

- ↵
- Brorsen KR,
- Yang Y,
- Pak MV,
- Hammes-Schiffer S

- ↵
- National Institute of Standard and Technology.

- ↵
- ↵
- ↵
- Kim M-C,
- Sim E,
- Burke K

- ↵
- Su NQ,
- Xu X

- ↵
- Su NQ,
- Xu X

- ↵
- ↵
- ↵
- Su NQ,
- Zhang IY,
- Xu X

- ↵
- Handy NC,
- Schaefer HF III

- ↵
- ↵
- Zhang IY,
- Rinke P,
- Scheffler M

_{2}/H_{2}^{+}challenge. New J Phys 18:073026. - ↵
- Zhang IY,
- Rinke P,
- Perdew JP,
- Scheffler M

- ↵
- Mezei PD,
- Csonka GI,
- Ruzsinszky A,
- Kállay M

- ↵
- ↵
- Langreth DC,
- Perdew JP

- ↵
- Gunnarsson O,
- Lundqvist BI

- ↵
- ↵
- Bader RFW

- ↵
- ↵
- ↵
- ↵
- Mezei PD,
- Csonka GI,
- Kállay M

- ↵
- Wang Y,
- Wang X,
- Truhlar DG,
- He X

^{6+}, and Ne^{8+}? J Chem Theory Comput 13:6068–6077. - ↵
- Mayer I,
- Pápai I,
- Bakó I,
- Nagy Á

- ↵
- Ryabinkin IG,
- Ospadov E,
- Staroverov VN

- ↵
- Goerigk L, et al.

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Physical Sciences