Poisoning density functional theory with benchmark sets of difficult systems †
Large benchmark sets like GMTKN55 [Goerigk et al., Phys. Chem. Chem. Phys., 2017, 19, 32184] let us analyse the performance of density functional theory over a diverse range of systems and bonding types. However, assessing over a large and diverse set can miss cases where approaches fail badly, and can give a misleading sense of security. To this end we introduce a series of 'poison' benchmark sets, P30-5, P30-10 and P30-20, comprising systems with up to 5, 10 and 20 atoms, respectively. These sets represent the most difficult-to-model systems in GMTKN55. We expect them to be useful in developing new approximations, identifying weak points in existing ones, and to aid in selecting appropriate DFAs for computational studies involving difficult physics, e.g. catalysis.It is becoming rare to find a physical chemistry paper that does not contain results from a density functional theory (DFT) calculation. DFT provides quantum chemical insights at low cost, by using density functional approximations (DFAs) to capture difficult quantum mechanics. [note, many authors erroneously use DFT to refer to DFAs] Popularity has brought with it diversity -there is now a 'zoo' of hundreds, if not thousands of DFAs to choose from. [1][2][3][4] The zoo continues to grow, as density functional developers produce new DFAs to solve outstanding problems and to cook new DFAs using new ingredients.
<!>How well does performance on WTMAD-2 predict performance on outliers?<!>Can we better understand density functional approximations by using outliers?
<p>Selecting the best DFA for a given study has thus become an onerous task, which has given rise to benchmarking studies that seek to make the task less onerous. Benchmarking helps to identify useful DFAs by scrutinizing their performance on relevant reactions. In broad terms, the process of benchmarking first requires computing,</p><p>for DFA d on 'reaction' I, which might be a traditional reaction or another important energy difference like an ionisation potential.</p><p>∆E d I is the energy difference calculated using the DFA. ∆E ref I is a reference value for the reaction computed using a high-level theory or obtained from appropriately modified experimental data.</p><p>a Qld Micro-and Nanotechnology Centre, Griffith University, Nathan, Qld 4111, Australia; E-mail: t.gould@griffith.edu.au Their absolute difference is thus a reasonable metric for the error of the given DFA on the given reaction.</p><p>The "general main group thermochemistry, kinetics, and noncovalent interactions" (GMTKN55) benchmark database 5 , consisting of ∼ 1500 reactions and ∼ 2500 individual energies, goes one step further than traditional benchmarking approaches. It condenses performance on a large benchmark set into a single number (WTMAD-2) for a given DFA, which is designed to assess a variety of predictions with different energy scales. GMTKN55 may thus be used to assess performance of existing DFAs, and also aid in optimising and scrutinising new approaches. Other sets, like MGCDB84 6 and MB16-43, 7 serve similar roles.</p><p>The WTMAD-2 scheme (see Sec 4 of GMTKN55 5 for details) defines an overall quality metric, WTMAD(d) := 1</p><p>using eq. ( 1) and weights, W I , that depend only on the benchmark set B(I) containing I. These weights normalise the deviations to ensure that sets with large energies, yet small relative errors, do not dominate over sets with small energies, yet relatively large and important errors. Thus, WTMAD-2 can be used to assess the effectiveness of a DFA, in practice.</p><p>However, a problem with any ranking protocol defined on a very large number of systems, such as WTMAD-2, is that the single number can mask systematic deficiencies of DFAs. For example, a DFA that very accurately predicts 1450 reactions can fail rather badly on the remaining 50 outliers without much statistical impact. One of the authors 8 previously showed that WTMAD-2 statistics could be reproduced imperfectly by just 50 carefully chosen reactions, and nearly perfectly by 150 reactions, using 'Diet-GMTKN55' benchmark sets. The ability to reduce the benchmark set size without impacting statistics might suggest that outliers are unlikely to be disastrous, since they are likely to have a greater impact on the smaller set. A more troubling possibility is that outliers simply did not make it into the 'diet' sets due to their low statistical importance.</p><p>When assessing DFAs, it would be useful to know how well they work for difficult cases, to test their overall robustness and reliability. The present work thus seeks to answer the following questions:</p><p>1. What reactions in GMTKN55 are most difficult for DFA to reproduce?</p><p>2. How well does performance on WTMAD-2 predict performance on outliers?</p><p>3. Can we better understand density functional approximations by using outliers?</p><p>Answering these questions will be the subject of the rest of this manuscript.</p><p>What reactions in GMTKN55 are most difficult for DFA to reproduce?</p><p>Answering the first question requires identifying what makes a reaction difficult. Any given DFA will struggle with some reactions -statistics basically guarantee it. A difficult reaction is therefore one which is difficult for a large number of DFAs to reproduce and which therefore poisons the statistics for a large number of DFAs. One therefore seeks systems that are outliers across many DFAs.</p><p>Identifying these reactions involves turning the usual benchmarking problem on its head. Eq. ( 2) assigns each given DFA, d, a (weighted) average of its performance across all reactions, I, which can be used to rank it -lower values indicated better quality. By contrast, the expression,</p><p>yields higher value for any reaction that is difficult for all N DFA DFAs to reproduce, and which thus poisons the most DFAs. The reactions may then be sorted from worst (largest Poison(I) value) to best (smallest Poison(I) value).</p><p>Selecting reactions with only the largest values of Poison(I) therefore defines a new benchmark set that is, by design, difficult for DFAs to reproduce. The smaller this set, the more difficult it will be for DFAs to reproduce the systems (by contrast, the entire set simply reproduces GMTKN55). However, smaller sets expose a potential deficiency in eq. ( 3) -it is based on absolute energies only and so may be biased to large systems with more bonds, since each bond can introduce error when computed using a DFA.</p><p>Maximal discriminatory power is therefore obtained by selecting a sufficiently small set of poison species, and ensuring it is not biased toward large systems. Poison30-N atom (P30-N) benchmark sets achieve this task, by using the thirty most difficult systems, as ranked by eq. ( 3), whose reactions contain at most N atom atoms. Thirty reactions (2% of the total GMTKN55 set) is chosen as a manageable number, with enough reactions for variety but few enough reactions to contain only outliers. The issue of bonds is avoided by using only reactions with up to a given number of atoms in its largest molecule. The present work settles on P30-5, P30-10 and P30-20, with five, ten and 20 atoms respectively. [One may also weight eq. ( 3) by W I before sorting, but defining a useful weight is not a trivial task.]</p><p>Figure 1 shows the molecules (but not atoms) used in P30-5. The other two sets are shown as Supplementary Figure 1 in the Supporting Information. Reactions in all poison sets come from just eight of the 55 sets in GMTKN55: (1) ALK8 5 ; (2) ALKBDE10 9 ; (3) BH76 [10][11][12] ; (4) DC13 [12][13][14][15][16][17][18][19][20][21][22][23] ; (5) G2RC 12,24 ; (6) MB16-43 7 ; (7) SIE4x4 5 ; and (8) W4-11 25 , with MB16-43, SIE4x4 and W4-11 making up the majority of the P30 sets in which they feature. Text descriptors in Figure 1 and its supporting counterparts show the original database, using these labels, and descriptions of each molecule from each original database.</p><p>The reactions identified for each P30 set can be rationalised by a combination of chemical and density-functional theorist intuition. The requirements to be selected in a poison set is to possess the highest error below a certain system size. The highest errors will therefore come from cases where (a) DFT performs particularly poorly, and/or (b) bond energies (and by extension errors) are particularly high, preferably both.</p><p>The three benchmark sets that feature prominently in the P30 sets match this intuition. Self-interaction error [26][27][28][29][30][31][32][33][34] is a well known DFA error which is often tested using dissociation of homodiatomic charged molecules. The SIE4x4 set is specifically designed to test for self-interaction error, so clearly satisfies (a), with charge imbalances common in these systems driving up energy differences satisfying (b) also. W4-11 is a set of total atomisation energies which satisfy requirement (b). Taking a closer look at some of the specific systems reveals a preference for homoneuclear molecules or multiple double bonds, systems known to exacerbate static/strong correlation [35][36][37][38] which is poorly described by DFAs and thus satisfy (a). MB16-43 is a set of randomly generated 'artificial molecules' and in not following the "narrow structural space of chemical intuition" 7 so it is not hard to image satisfying requirements (a) and (b). The remaining contributors to the P30 sets from other benchmark sets also fit into the arguments given here.</p><!><p>With the P30-N sets established, the second question may now be addressed. First, the unweighted MAD,</p><p>gives an error metic for each DFA, d, on poison set P30-N. Then, the error of each DFA on P30-N may be compared against the WTMAD-2 value for the same DFA using eq. ( 2). Following GMTKN55, DFAs are divided into different rungs of "Jacob's ladder". 39 The lowest rung is GGAs, which are function-als of the density, n and its gradient, |∇n|, only. The next rung is meta-GGAs (MGGAs) which introduce the kinetic energy density, τ, and Laplacians of the density, ∇ 2 n. The next rung is hybrids, which introduce some Hartree-Fock (HF) exchange into the functional -we do not distinguish between range-separated and traditional hybrids. The highest rung is double hybrids, (DHDF) which add a wave-function based correlation energy (usually second order Møller Plesset perturbation theory) into the mixture. As a rule of thumb, higher rungs are expected to be more accurate than lower ones, but also take longer to compute. A D3(BJ) dispersion correction is used in almost all cases due to its ability to improve most DFAs in GMTKN55, with exceptions being approaches where D3(0) or another dispersion correction are more appropriate. Choices are made consistent with Sec 3.2 of Sec 4 of GMTKN55 5 . Only the ten best methods of each rung are reported, when there are more than ten available. The reported hybrids vary significantly with the number of atoms in the poison set. However, this most likely reflects the fact that there are 47 hybrids to choose from, many with similar performance, versus 17 GGAs in the DFA suite. All 7 double hybrids and 8 meta-GGAs are reported. Supplementary Figure 2 shows errors for all DFAs broken down by set and rung.</p><p>To go beyond visual comparisons, a Kendall rank correlation coefficient 40 is computed for each of the rungs (including DFAs not shown), using, Four main conclusions may be drawn from the Kendall coefficients and Figure 2: i) hybrids and GGAs have minimal correlation between the poison sets and WTMAD-2; ii) double hybrids are reliable in the sense that they performance on difficult (poison) systems reflects their overall performance, as indicated by relative large values of τ; iii) meta-GGAs tend to sacrifice quality on poison reactions for better overall accuracy, as indicated by negative values of τ; iv) correlation is much lower in P30-20 than the other two sets, as indicated by |τ| approaching zero.</p><!><p>Finally, we address the question of what we can learn from the poison sets. Many of the most accurate DFAs on GMTKN55 involve multiple empirical numbers of parameters that are optimized on large benchmark sets. Unsurprisingly, these also perform well on the poison sets, albeit less well and less consistently than on the complete set. Unfortunately, the large number of pa-rameters make it difficult to understand sources of errors.</p><p>PBE0 41,42 and B3LYP 43,44 offer a simpler level of empiricism, which do allow analysis. PBE0 may be considered to have one empirical parameter being the fraction of Hartree-Fock exchange used in the hybrid,</p><p>and which is set to 0.25 in regular PBE0. B3LYP enhances this with a second parameter for the mixture of LDA and GGA and uses α ≈ 0.2. Since PBE0 has just one parameter we investigate it further. Barring exceptional cases, each reaction will have a value, α 0 , such that the errors of HF and PBE cancel out. That is, we can find α 0 (I) for each reaction I such that Err(PBE0 α 0 (I) , I) = 0.</p><p>We compute errors as a function of α for all reactions in P30-5. This lets us determine α 0 , which we plot in the left panel of Figure 3. The right plot shows errors as a function of α, broken down by subset (i.e. MAD within each subset) and over the full poison set -we restrict to subsets with more than two reactions present in P30-5. In some cases, calculations with 0 ≤ α ≤ 1 did not contain the optimal value, so we extrapolated α 0 -see Supporting Information. G2RC:11 (exchange of Si and C in SiF 4 and CCl 4 ) is particularly bad, with a predicted and highly unphysical We return to the prominent benchmark sets SIE4x4 and W4-11 in order to discuss patterns in α 0 revealed by Figure 3. SIE4x4 consistently gives α 0 ≈ 0.7. Self-interaction errors are not present in HF theory and so these larger α values effectively offset errors seen in SIE4x4. W4-11 is an atomisation benchmark set yielding α 0 scattered around 0.25. This is consistent with the ideal α observed in most hybrid studies which focus on atomization energies as the benchmark, especially the earliest hybrid works. 42,43,45 Remarkably, the optimal α for all full poison sets, that is P30-5, P30-10 and P30-20, is also close to the PBE0 value of 0.25. This can be seen from the right plot of Figure 3 and Supplementary Figure 3. However, we also note that the optimal value of α is clearly an average, and that α 0 is very system specific. Indeed, much has been made of the seemingly fortuitous one quarter balance of exact exchange, but with little consensus emerging in the last two decades. Optimally tuned hybrid functionals,  which identify an optimal α for a given reaction using accessible criteria, are an increasingly common way to avoid the 'scatter' in α 0 . This analysis may offer insights into the good behaviour of the recently introduced DM21 DFA, 69 which is based on a machinelearned local-hybrid ansatz in which α varies in space. Because α can vary in space, DM21 is able to 'see' local conditions and adjust the weight accordingly. This feature cannot be replicated by global hybrids where α is a unique constant -all hybrid DFA reported except DM21 are global in nature. As is revealed in Table 1, DM21 out-performs global hybrids (by a factor of five on P30-20, versus ωB97x-V 70 ) for difficult 'poison' cases, without sacrificing its strong overall performance (here illustrated by WTMAD-2) which is as good as top-performing (on WTMAD-2) ωB97x-V and better than second-ranked M052X.</p><p>Interestingly, PBE0 (one parameter) performs much better on the largest poison set than the highly-empirical and accurate ωB97x-V (around 15 parameters) and M052X 71 (around 20 parameters), and the popular B3LYP (two parameters). [Note, DM21 has effectively thousands of parameters] The reason for this unusual success should be investigated thoroughly as it may offer strategies to improve robustness of global hybrid DFAs.</p><p>To conclude, this work introduced several new benchmark sets of difficult 'poison' reactions. The three sets, P30-5, P30-10 and P30-20, consist of thirty reactions involving systems with up to five, ten and 20 atoms, respectively. Each set is composed of the reactions in GMTKN55 that are most difficult for DFAs to reproduce, and which likely offer the most challenging test for any DFA. They thus provide benchmark sets that are particularly difficult to model, and which serve as a stringent test of any DFA.</p><p>With the exception of double hybrid DFAs, which are quite expensive to compute, a high overall ranking on WTMAD-2 did not necessarily mean good success on P30-N. This means that overall success is only weakly correlated with success on difficult systems, especially as they become larger and especially for meta-GGAs, which should therefore be used with caution. Selection of a DFA using WTMAD-2 (or, we expect, similar metrics) may therefore be of little predictive value when studying difficult reaction processes, such as those involved in complex catalysis.</p>