Technical CommentsTheoretical Chemistry

Response to Comment on “Density functional theory is straying from the path toward the exact functional”

See allHide authors and affiliations

Science  05 May 2017:
Vol. 356, Issue 6337, pp. 496
DOI: 10.1126/science.aam9550

Abstract

Kepp argues in his Comment, among other concerns, that the atomic densities we have considered are not relevant to molecular bonding. However, this does not change the main conclusion of our study, that unconstrained fitting of flexible functional forms can make a density functional more interpolative but less widely predictive.

Kepp argues in his Comment (1) that our Report (2) does not properly assess the accuracy of approximate functionals for molecules. Many studies have been published to assess the accuracy of such functionals for various applications. We have, instead, investigated the ability of these functionals to make predictions for electron densities—the property employed in computations of all other properties but almost never used to fit or assess functionals. We deliberately chose atomic systems because they are somewhat different from molecules and because they should be well described by the tested functional forms. We found that nonempirical or few-parameter empirical functionals were much more predictive than flexible (highly unconstrained) multiparameter empirical functionals. We will respond in detail to the Comment after clarifying the methodological issue in density functional theory (DFT) and more generally.

There is a systematic, nonempirical way to improve approximations to the exact density functional for the exchange-correlation energy (3, 4): (i) Prove the existence of the exact functional and derive exact formal expressions for it. (ii) Discover mathematical properties of the exact functional (exact constraints), which include limits, scaling relations, equalities, and bounds. (iii) Develop approximate but computationally tractable forms for the approximations at various levels of flexibility, including the local spin density approximation (LSDA), generalized gradient approximations (GGAs), meta-GGAs, and hybrids tested in our paper. (iv) Impose the “exact constraints” from step ii on each form, as appropriate. (v) If a form still retains some flexibility, fit it to energies/densities of appropriate norms (4), systems for which the form can be expected to be highly accurate. (Bonded systems are not appropriate norms, because in them the standard approximations display an understood but uncontrollable error cancellation between exchange and correlation. Fitting to bonded systems could make a functional more interpolative but less widely predictive, in the sense of the next paragraph.) The period 1965 to 1979 focused on step i, the 1980s on step ii, 1965 to 1998 on step iii, and the 1970s to the present on step iv, whereas the first appropriate norm for step v was the uniform electron gas (1965). A recent example of steps i to v is the SCAN (strongly constrained and appropriately normed) meta-GGA (4), which respects all 17 exact constraints that a meta-GGA can and works well for diversely bonded molecules and materials [e.g., reference 5 in (2)]. Because SCAN was fitted only to energies of nonbonded systems, every bonding description from it is a genuine prediction. On the meta-GGA level, SCAN performs slightly better than the multiparameter empirical M06-L (5) on a large database of main-group chemical energies (with weighted mean absolute deviations of 4.1 and 4.9 kcal/mol, respectively) (6), and we found it to be much better for the electron density (2).

A shortcut that has become increasingly popular is to skip or downplay steps ii, iv, and v and to use the bonded systems, which are often of greatest interest to the user, for fitting. There have even been suggestions, as in (7), to drop exact constraints previously incorporated in the functional forms to get better fits. In an era of big data and machine learning (8), it is tempting to think that one only needs to expand the functional flexibility and the fitting sets. Indeed, that may be enough to interpolate between similar systems and properties, but our results (2) suggest that only incorporation of the exact constraints can result in reliable prediction over the enormous space of possible molecules, materials, and properties. The interpolation/prediction dichotomy is of course greater for unconstrained machine-learned functionals than it is for the 128 functionals that we studied, which all at least start from physical insights.

An important but often implicit exact constraint on an LSDA, GGA, or meta-GGA is the smoothness of its enhancement factor over exchange-only spin-unpolarized LSDA. Although there is no proof that the enhancement factor must be smooth, Occam’s razor demands that one must first assume and test that simpler hypothesis. Many functionals that yielded the best densities and density derivatives in our tests are smooth, whereas flexible empirical functionals that yielded the worst errors display oscillations in the enhancement factor [see figure 1 in (9)] that could produce oscillations of the exchange-correlation potential and thus of the density [figures S1 to S3 in (2)]. Oscillations of the density are magnified in its derivatives. Fitting to densities and their derivatives could also smooth the enhancement factor of an approximate functional.

Would the oscillations in the densities from the multiparameter empirical functionals be removed by using a finer grid? Our calculations were performed (2) using the Ultrafine or the JANS=2 grids, with 99 and 155 radial points, respectively. We have now tried a grid with 999 radial and 974 angular points, finding no visible change for most functionals, and visible but insignificant changes (less than 0.1 units of median-normalized absolute error) for M06 and M11-L, on the systems Ne, Be, Ne6+, and F5+ from (2).

After this preamble, we respond to the first two sentences of the Comment: While the GGA potentials described by Cruz et al. made larger relative errors than the GGA total energies, the step from LSDA to GGA still improved both densities (2) and energies. Now, we respond to each of the numbered points in (1).

1) Our conclusion is that the increasingly used methodology of fitting flexible functionals without full regard for exact constraints strays from the path toward the exact functional. Although in our study, many of these functionals are from the Minnesota research group, the trend is more general. We tested a very large, but not exhaustive, set of functionals. More will be included in future work when we are able to access them.

2) Our “path toward the exact functional” is a methodology and thus not plottable. It is nevertheless true that properly constrained functionals improve in both energy and electron density from LSDAs to GGAs to meta-GGAs to hybrids, whereas unconstrained functionals can worsen the density in this progression. A direct comparison of errors in energies and electron densities is not useful, however, because (i) along with the density error, which we study in our article, there is always a functional error (an error made by the functional on the exact electron density), which interferes inter se (10); and (ii) energies (except the total energies) are relative values, which correspond to the difference between two quantum systems, each of which has its own density error and functional error. Thus, good energies do not necessarily mean good densities, and this was one of the main points of our article: Good performance for energies is never enough for a good functional because it can arise from an uncontrollable error cancellation.

The Comment raises a question about the influence of errors in our electron density descriptors on the energy. All three of our density descriptors are relevant to the energy, but we cannot at present answer more quantitatively. Whatever the answer, we have shown that features of the density that are predicted well by constrained functionals (and would be exact for the exact one) are seriously worsened by fitting a flexible functional form to molecular energies.

3) We agree that Hartree-Fock performs well for our two-electron ions (with very little correlation) and not for our four-electron systems. Its ranking in our study is only mediocre. Our top performers—few-parameter hybrid functionals—are on top not because they are very accurate for some systems but because they are not too inaccurate for any of them. This was the idea behind sorting of the functionals by their maximum absolute errors.

4) We do not believe in a universal 25% of exact exchange in a hybrid functional. A model (11) for the coupling-constant dependence of the atomization energy suggests that 25% is about right for a correction to the PBE GGA when applied to a moderate-energy-gap molecule or solid. Figure S5 in (2) shows that the large-gap Ne atom requires about 50% of exact exchange to optimally correct PBE. Although the optimum percentage is system-specific, we must choose it to be a constant or the method will not be size-consistent. A 25% correction to PBE, as in the hybrid PBE0, is not bad for all our studied systems.

5) Our atoms and atomic ions were not intended as proxies for molecules. Instead, they are appropriate norms (4): easy cases for the considered mostly semilocal functionals. Functionals that make unnecessarily large errors for the densities of these systems have internal problems that can be reflected in applications to more challenging systems, such as molecules. As an illustration, we have calculated the electron densities of a water molecule (at the experimental gas-phase geometry taken from the National Institute of Standards and Technology database) with two density functionals, M11 (which gave less accurate densities for atoms) and PBE0 (which gave more accurate ones) and then subtracted from them the nearly exact CCSD-full density (Fig. 1).

Fig. 1 Comparison of the errors in the density of the water molecule between M11 and PBE0.

(Left) ρM11 – ρCCSD-full. (Right) ρPBE0 – ρCCSD-full. Isosurfaces are drawn at Δρ = 0.01 e/bohr3 (blue) and Δρ = –0.01 e/bohr3 (red). The basis set was aug-pc-3 (15). Near the oxygen lone pairs, the errors reach −0.0151 e/bohr3 (−1.5%) for M11 and −0.0092 e/bohr3 (−0.9%) for PBE0.

The Comment suggests that the large errors we found for the density of the neon atom may arise from relativistic s-shell contraction, but that cannot be so because our approximate and reference densities are both nonrelativistic.

In figure 1 of (1), the density and energy errors are roughly proportional to one another because each scales like the nuclear charge at fixed electron number for a given functional. That figure also shows that the multiparameter empirical hybrid M06-2X functional is excellent for the two-electron removal energies of our four-electron ions and also makes a remarkably small median-normalized absolute error in the density. We agree that, if one neglects its errors in the derivatives of the density, M06-2X would be one of the best-performing (instead of one of the worst-performing) functionals in our study.

Overall, Minnesota functionals deserve major credit for their important contributions to DFT development. M06-L was the first functional to show (5) that a meta-GGA can describe intermediate-range van der Waals interactions. More generally, Minnesota functionals have shown that carefully fitted flexible forms can make accurate predictions for systems and properties similar to those on which they were fitted [e.g., for diverse possible transition states of the SpnF-catalyzed formal Diels-Alder cycloaddition (12)]. However, their neglect of some exact constraints can lead to uncontrolled errors for systems and properties that are very different [e.g., in (2), and for Au bonding, as recently noted by Kepp (13)]. [See (8, 14) for excellent discussions of such fitting issues.] Concerning the performance of M06-2X, it should be noted that, as its authors carefully stated in their abstract (5), M06-2X (a hybrid with double the exact exchange of M06) is not a general-purpose functional but one designed specifically for main-group molecules. In that limited but important domain, it may be highly accurate for energies, but it is off the path toward the exact universal functional. Closeness to the exact functional requires not only accurate energies on broad data sets but also accurate densities and satisfaction of exact constraints.

References and Notes

Acknowledgments: M.G.M., I.S.B., and K.A.L. are grateful to Russian Science Foundation grant 14-13-00884 for financial support. J.P.P. and J.S. were supported by the U.S. National Science Foundation under grant DMR-1607868. M.G.M. acknowledges the computational resources provided by the Supercomputing Center of Lomonosov Moscow State University: the Lomonosov supercomputer.
View Abstract

Navigate This Article