## Abstract

Medvedev *et al*. (Reports, 6 January 2017, p. 49) argue that recent density functionals stray from the path toward exactness. This conclusion rests on very compact 1s^{2} and 1s^{2}2s^{2} systems favored by the Hartree-Fock picture. Comparison to actual energies for the same systems indicates that the “straying” is not chemically relevant and is at best specific to the studied dense systems.

Medvedev *et al*. (*1*) point out that electron densities ρ and energies *E*[ρ] computed with density functional theory (DFT) do not always increase in accuracy together. Cruz *et al*. (*2*) stated the problem in 1998 as “functionals which yield highly accurate energies often produce potentials which differ markedly from the exact ones.” Medvedev *et al*. put errors in ρ on a time scale and show a trend of improvement impaired by nine specific recent functionals with reported high accuracy of *E*[ρ] for diverse systems. The inverse relationship in Fig. 1B may suggest an overfitting problem on the path toward universality, where both ρ and *E*[ρ] should become increasingly accurate; off this track, accurate energies with inaccurate densities would seem successful only until applied outside the parameterization range. Some comments seem warranted.

1) Of the nine specific functionals that deviate from the “path,” almost all are from 2011–2012 and all are from one specific research group; other recent functionals perform well in the trend, and the only two functionals from 2014–2015 are on-path. With two functionals from 2015, none from 2014, and three from 2013, the recent history seems undersampled; various post-2011 functionals by other groups have not been included (*3*–*7*). Thus, whereas some recent functionals from one research group have apparently sacrificed some accuracy in 1s^{2} and 1s^{2}2s^{2} systems for accuracy in diverse molecular energies, arguing that DFT deviates from the path seems an overgeneralization.

2) Another concern is whether the functionals are actually on a “path,” as no direct comparison of *E*[ρ] and ρ was done; the errors in *E* were from general benchmarks of diverse molecules (*8*). There is only a path if the errors of both *E*[ρ] and ρ decrease together for the same systems, and the studied systems are very distinct. Also, the authors used maximum errors after normalization of both ρ, its gradient, and Laplacian for ranking, which gives specific weights to terms that perhaps do not reflect their importance, as would be measured by their relative effect on *E*.

3) The Hartree-Fock (HF) method has a root mean square deviation (RMSD) for ρ of only 0.049 for the 1s^{2} systems (B^{3+}, C^{4+}, N^{5+}, O^{6+}, F^{7+}, and Ne^{8+}) [data S3 of (*1*)]. Six of the 14 systems studied (43%) are of this type. The high accuracy of HF is specific to systems with 2N^{2} valence electrons (the octet rule), where N is the period number [e.g., figure 2 in (*1*), where Ne requires much more HF exchange, as do the 1s^{2} systems]. If one leaves out the six two-electron systems, HF exhibits worse performance than many functionals (average RMSD of ρ = 1.81 without two-electron systems, 0.92 with). Thus, the choice of benchmark systems favors the HF picture and is not a reasonable choice of norm. Accordingly, the top performers are all hybrid functionals and deviations from “exactness” are at best specific to these systems.

4) Similarly, figure 2 in (*1*) reports the maximum error of ρ, its gradient, and Laplacian (the inclusion of the latter affects the ranking, e.g., of M06-2X); inclusion of the six 1s^{2} systems would reveal the high HF demands of the 1s^{2} configurations, and for other portions of the periodic table, smaller HF percentages are required (*9*, *10*), yet in the 1s^{2}2s^{2} systems the gap between virtual and occupied orbitals justifies 25%. Thus, a version of figure 2 in (*1*) with all systems included would indicate that the role of the 25% is diminished substantially.

5) Thirteen of the systems have 1s^{2} or 1s^{2}2s^{2} configuration, and 10 of the 14 studied ions have a charge ranging from +3 to +8, representing extremely compact ρ with large dynamic (but no static) correlation—namely, the large improvement by MP4 over MP2 [data S1 of (*1*)]. Such compact densities are not found in ordinary reaction chemistry because they require tens of eV to generate; it is thus questionable whether this very compact ρ regime is chemically relevant.

To address points 1 to 5, because energy is a state function, the quality of *E*[ρ] can be probed by comparing with ionization potentials (IP) from the National Institute of Standards and Technology database—e.g., *E*[ρ] of B^{3+} and B^{+} can be probed by the second and third experimental IP of boron (dication energies cancel out). *E*(B^{3+}) − *E*(B^{+}) = IP3(B) + IP2(B) = 37.931 eV + 25.155 eV = 63.085 eV(1)This experimental energy corresponds to removal of both 2s electrons from the 1s^{2}2s^{2} configurations, with a trend of increasing charge. Comparing with *E*[ρ] directly reveals whether errors in ρ have chemical relevance and whether there is a relationship between errors *E*[ρ] and ρ implying a “path” toward universality and, accordingly, a deviation from such path, as claimed.

Computations were carried out with Turbomole 7.0 (*11*) for *E*(B^{3+}) − *E*(B^{+}) = IP3(B) + IP2(B) = 63.085 eV; *E*(C^{4+}) − *E*(C^{2+}) = IP4(C) + IP3(C) = 112.381 eV; *E*(N^{5+}) − *E*(N^{3+}) = IP5(N) + IP4(N) = 175.364 eV; *E*(O^{6+}) − *E*(O^{4+}) = IP6(O) + IP5(O) = 252.018 eV; *E*(F^{7+}) − *E*(F^{5+}) = IP7(F) + IP6(F) = 342.350 eV; and *E*(Ne^{8+}) − *E*(Ne^{6+}) = IP8(Ne) + IP7(Ne) = 446.368 eV. This sampling covers 12 of the 14 systems using the same aug-cc-pwCV5Z basis set, tight densities, energies, and grids. For illustration, PBE0, TPSSh and TPSS, B3LYP, BHLYP, BP86, M06 and SVWN, M06-2X, HF, MP2, and CCSD were studied. They spread across the ranking by Medvedev *et al*. CCSD(T) was also included because CCSD, although a full-CI method and thus exact nonrelativistically for 1s^{2} systems, may miss some core-valence correlation of the 4- and 10-electron systems.

Figure 1A (nonrelativistic) and 1B (corrected for relativistic effects) show that relativistic effects grow with charge, as 1s-electrons are accelerated. Relativistic stabilization and contraction of the s-shells favor the 1s^{2}2s^{2} systems over 1s^{2} systems. Due to zero spin and angular momentum, scalar relativistic corrections recover this effect (Fig. 1B) and are >0.6 eV for the neon systems [the neon systems have the largest errors in (*1*), probably because ρ is relativistically contracted]. Relativistic corrected CCSD(T) and CCSD energies are within 0.03 eV (~3 kJ/mol) of experiment. Accordingly, the exact density functional methodology would provide exact energies to within 3 kJ/mol if applied with this basis set and relativistic correction. Thus, we can compare the density functionals now also in the energy regime, *E*[ρ].

HF errors in energy exceed 3 eV for neon systems (Fig. 1B). Local functionals M06 and SVWN produce errors almost as large as HF. Most other functionals perform similarly, although B3LYP and M06-2X perform distinctly better. Only the first bar represents chemical relevance, as net atomic charges in molecules rarely exceed 2 under ambient conditions even for highly charged molecules. For chemically relevant boron, all DFT methods perform better than MP2, which only becomes more accurate in the very compact, highly charged ions. The error of B3LYP is 0.03 eV, and the worst-performing functionals (PBE0 and BP86) show 0.27 to 0.28 eV. Thus, the extremely compact regime mostly studied by Medvedev *et al*. is probably not chemically relevant, yet clearly affects the ranking.

To produce consistent paths toward exactness, one has to study *E* and ρ for the same systems. To this end, the RMSD of ρ of the largest 1s^{2}2s^{2} ions from Medvedev *et al*. [data S4 in (*1*)] was compared to errors in the energy of removing the two 2s^{2} electrons. Figure 1C reveals almost perfectly linear relationships. Because all the energies are for iso-electronic conversions, this relation reflects a sensitivity to charge, which accelerates the electrons and increases the kinetic energy and correlation in the very compact systems. Most DFT methods and MP2 follow a “path” of accuracy with errors in energy growing with errors in ρ (coefficients of −0.77 to −0.86). M06-2X errors in *E*[ρ] increase slowly with ρ, whereas the local SVWN and HF energies deteriorate much more rapidly as ρ becomes compact. Notably, M06-2X is very “exact” when put on actual *E*[ρ] paths and much more exact than MP2, PBE0, or TPSSh. In table 2 in (*1*), M06-2X was ranked very low mainly because of the Laplacian of ρ and thus claimed to be off-path, despite *E*[ρ] and ρ being excellently on-path (Fig. 1C).

In conclusion, the poor performance of some recent functionals for very compact densities of highly charged closed-shell systems does not imply that they are less exact. Because hybrid functionals are favored by system choice, they perform best in the test. More interestingly, functionals show distinct error relationships between ρ and *E*[ρ] (Fig. 1C), with exactness being represented by CCSD(T) in the lower right corner. These relationships are on actual paths and are likely to be focus points in targeting exact functionals, but they would have to pertain to densities that are more chemically relevant and diverse.

## Supplementary Materials

www.sciencemag.org/cgi/content/full/356/6337/496-b/DC1

Data S1 and S2

## References and Notes

**Acknowledgments:**The author acknowledges computer time at the Centre for Scientific Computing at Aarhus University.