The Elephant in the Room of Density Functional Theory Calculations
Abstract
Using multiwavelets, we have obtained total energies and corresponding atomization energies for the GGAPBE and hybridPBE0 density functionals for a test set of 211 molecules with an unprecedented and guaranteed Hartree accuracy. These quasiexact references allow us to quantify the accuracy of standard allelectron basis sets that are believed to be highly accurate for molecules, such as Gaussiantype orbitals (GTOs), allelectron numeric atomcentered orbitals (NAOs) and fullpotential augmented plane wave (APW) methods. We show that NAOs are able to achieve the socalled chemical accuracy (1 kcal/mol) for the typical basis set sizes used in applications, for both total and atomization energies. For GTOs, a triplezeta quality basis has mean errors of 10 kcal/mol in total energies, while chemical accuracy is almost reached for a quintuplezeta basis. Due to systematic error cancellations, atomization energy errors are reduced by almost an order of magnitude, placing chemical accuracy within reach also for medium to large GTO bases, albeit with significant outliers. In order to check the accuracy of the computed densities, we have also investigated the dipole moments, where in general, only the largest NAO and GTO bases are able to yield errors below 0.01 Debye. The observed errors are similar across the different functionals considered here.
Electronic structure calculations are nowadays employed by a large and steadily growing community, spanning condensed matter physics, physical chemistry, material science, biochemistry and molecular biology, geophysics and astrophysics. Such a popularity is in large part due to the development of Density Functional Theory (DFT) methods,Hohenberg and Kohn (1964) in their Kohn–Sham (KS) formulation.Kohn and Sham (1965)
Although the exact energy functional of DFT is unknown, many approximate functionals offer an excellent compromise between accuracy and numerical cost, rivaling often the accuracy that can be obtained with correlated methods, such as CoupledCluster Singles Doubles (CCSD). Goerigk and Grimme (2010); Karton et al. (2008); Zhao and Truhlar (2005) During the last decades, extensive efforts have been undertaken to provide ever more accurate approximations to the exact ExchangeCorrelation (XC) functional. Perdew et al. (1996) This quest for higher accuracy is conceptually captured by John Perdew’s Jacob’s ladder analogy, Perdew and Schmidt (2001) leading to the heaven of chemical accuracy: errors of 1 kcal/mol or less in atomization energies and other energy differences that are of primary interest in chemistry and solid state physics. Rungs on this ladder are the Local Density Approximation (LDA), the Generalized Gradient Approximation (GGA), Li et al. (2015) metaGGAs, Sun et al. (2015) hybrid and double hybrid functionals. Su and Xu (2015) The best modern XC functionals come fairly close to this target, with errors of a few kcal/mol on a wide range of energetic properties relative to experiment, including atomic and molecular energies, bond energies, excitation and isomerization energies and reaction barriers, for maingroup elements as well as transition metals and solids. Peverati and Truhlar (2014); Mardirossian and HeadGordon (2016)
The closer we get to chemical accuracy, the more important becomes the identification of errors due to various other, algorithmic approximations – basis sets, integration grids and pseudopotentials, Singh and Nordstrom (2006); Goedecker and Maschke (1992); Willand et al. (2013) to cite a few – which can lead to comparable or even larger errors, but their influence is hard to quantify. The importance of this issue has recently been highlighted within the solid state community, with a substantial effort to assess the influence of such approximations on the accuracy. Lejaeghere et al. Lejaeghere et al. (2016) compared the GGAPBEPerdew et al. (1996) calculated equations of state for 71 elemental crystals from 15 different widely used DFT codes, employing both allelectron methods as well as 40 different pseudopotential sets. For the equation of state, most DFT codes agree within error bars that are comparable to those of experiment, irrespective of the basisset choice: allelectron numeric atomcentered orbitals, augmented plane wave (APW) methods or plane waves with pseudopotentials.
APW methods Singh and Nordstrom (2006) are widely believed to be highly accurate, but contain several parameters which are difficult to adjust and which can influence the results in a more or less erratic way. Hence, the magnitude of the error cannot be rigorously quantified without an external reference. Similar limitations exist for atomization energies of molecules obtained with Gaussiantype orbital (GTO) and NAO basis sets: both bases cannot be systematically enlarged to achieve completeness in the sense, and within standardized basis sets, the convergence to the exact result cannot be achieved. Additionally, for larger systems, linear dependency issues can limit the ability of these basis sets to achieve complete convergence.Moncrieff and Wilson (2005)
The basic mathematical formalism for KS DFT calculations leads to a selfconsistent threedimensional partial differential equation. What makes the solution of this equation so challenging are the accuracy requirements for the physically and chemically relevant energy differences. For instance, the atomization energy of the largest molecule (SiCl) in our data set is less than 1 Ha, but it is computed as a difference of energies in the order of 2000 Ha. Hence we need at least 7 correct decimal places in the total energy of the molecule to get the atomization energy within chemical accuracy of 1 kcal/mol milliHartree (mHa), and even 10 decimal places for microHartree (Ha) accuracy.
For isolated atoms, and using the appropriate numerical techniques, the associated manyparticle problem can be solved with essentially arbitrary numerical precision. Virtually converged LDA energies for spherical atoms are available in the NIST data base. Kotochigova et al. (1997) For a few dimers, highly accurate energies have been calculated, Ballone and Galli (1990) and an attempt to obtain total energies free of basis set error was made also for general molecular systems, Becke (1992) but the accuracy of this approach seems to be limited to around 1 kcal/mol. Similar accuracies were achieved for solids with semicardinal wavelets. Daykov et al. (2003)
In spite of all the progress that has been made in numerical techniques to harness the power of quantum mechanical theory and simulations, none of the traditional techniques is able to furnish, unambiguously, atomization energies for molecules with arbitrary numerical precision. A straightforward, uniform grid or Fourier transformbased approach is ruled out since it is impossible to provide sufficient resolution for the rapidly varying wave functions near the nucleus. Other basis set techniques are critically hampered by nonorthogonality, which leads to inevitable algebraic illconditioning problems at small but finite residual precision. Moncrieff and Wilson (2005) Because of these problems a large part of the community resorts to pseudopotentials methods, Singh and Nordstrom (2006) where the potential is replaced by a smoother potential that retains approximately the same physical properties of the allelectron atom. The smooth pseudopotential then allows to obtain arbitrarily high accuracy with systematic basis sets such as plane waves. The limitation is that the pseudopotential introduces an approximation error, the magnitude of which is hard to quantify.
The current de facto standard technique to assess errors of different methods does not rely on an absolute reference: errors are instead estimated by comparing results obtained with increasingly large bases. NIS (2015); Papas and Schaefer (2006) The development of Multiwavelet methods Harrison et al. (2003, 2004); Yanai et al. (2004a, b) have fundamentally changed the situation. Multiwavelets are systematic, adaptive, and can be employed in allelectron calculations. With this approach, it is now possible to achieve allelectron energies with arbitrarily small errors.
In the present work, we use MWs to obtain error bars of less than a Hartree (Ha) in the atomization energies for a large test set of 211 molecules with standard DFT functionals. We focus on three widely used and well established functionals, LDASVWN5, Vosko et al. (1980) GGAPBE, Perdew et al. (1996) as well as hybridPBE0. Adamo and Barone (1999) PBE and PBE0 are both relatively accurate for atomization energies, Paier et al. (2005) and have stood the test of time. Medvedev et al. (2017) Our MW results provide quasiexact reference values that can be employed to quantify the accuracy of standard basis sets, such as GTO, NAO, and APW methods, as well as of novel approaches based for instance on finite element methods Pask and Sterne (2005); Losilla and Sundholm (2012); Toivanen et al. (2015); Parkkinen et al. (2017) or discontinuous Galerkin methods. Lin et al. (2012)
Realspace methods have a long history in computational chemistry and have been used for benchmarking purposes for decades. Frediani and Sundholm (2015) However, because of the socalled curse of dimensionality, the naïve numerical treatment of molecular systems is prohibitively expensive, and its applicability relies on high symmetry to reduce the dimensionality. Sundholm and Olsen (1994) The multiscale nature of the problem renders the traditional uniform grid discretization highly inefficient, unless the problematic nuclear region is treated separately, e.g. by means of pseudopotentials. The mathematical theory to solve these issues was developed in the ’90s, when Alpert introduced the MW basis, allowing for nonuniform grids with strict control of the discretization error, as well as sparse representations of a range of physically important operators, with high and controllable precision. Alpert et al. (1993); Alpert (1999)
Alpert’s construction starts from a standard polynomial basis of order , such as the Legendre or the Interpolating polynomials, rescaled and orthonormalized on the unit interval . Then, an orthonormal scaling basis at refinement level is constructed by dilation and translation of the original basis functions , where is the th polynomial in the interval at scale . The set of scaling functions on all translations at scale defines the scaling space , and in this way a ladder of spaces is constructed such that the complete limit can be approached in a systematic manner:
(1) 
The wavelet spaces are simply the orthogonal complement of two subsequent scaling spaces and :
(2) 
Completeness in the sense can be achieved both by increasing the polynomial order (larger ) of the basis and by increasing the refinement in the ladder of spaces (larger ).
Two additional properties are essential in achieving fast and accurate algorithms: the vanishing moments of the wavelet functions and the disjoint support of the scaling and wavelet functions. The former leads to fast convergence in the representation of smooth functions and narrowbanded operators, whereas the latter enables simple algorithms for adaptive refinement of the underlying numerical grid, which is essential to limit storage requirements. The extension to several dimensions is achieved by standard tensorproduct methods; to minimize the impact of the curse of dimensionality, it is necessary to apply operators in a separated form, Beylkin and Mohlenkamp (2002); Beylkin et al. (2007) by rewriting the full multidimensional operator as a product of onedimensional contributions. For many important operators such a separation is not exact, but Beylkin and coworkers have shown that it can be achieved to any predefined precision as an expansion in Gaussian functions. Beylkin and Mohlenkamp (2002); Fann et al. (2004); Beylkin (2005) To apply such operators in multiple dimensions, and simultaneously retain the local adaptivity in the representation of functions, it is essential to employ the nonstandard form of operators, Gines et al. (1998); Beylkin et al. (2008) which, in contrast to the standard one, allows to decouple length scales.
These combined efforts (MW representation of functions and operators, separable operator representations, nonstandard form of operators) made the accurate application of several important convolution operators efficient in three dimensions:
(3) 
Among such operators are the Poisson () and the BoundState Helmholtz (BSH) () kernels:
(4) 
This mathematical framework was introduced to the computational chemistry community in the mid 2000’s by Harrison and coworkers. Harrison et al. (2003, 2004); Yanai et al. (2004a, b) They demonstrated that MWs could be employed to solve the KS equations in their integral reformulation: Kalos (1962)
(5) 
where the ’s are the KS orbitals and the potential operator includes external (nuclear), Hartree, exchange and correlation contributions, while the kinetic operator and the orbital energies are included in the BSH operator as , with . The ordinary KS equations , where is the KS Hamiltonian, can be recovered by recalling that . Such an integral formalism, when combined with a Krylov subspace accelerator, Harrison (2004) leads to fast and robust convergence of the fixpoint iteration of Eq. (5).
The integral formulation, in combination with MWs, provides unprecedented accuracy for allelectron calculations, without relying on molecular symmetry, Yanai et al. (2004b, a); Harrison et al. (2004) and has also been extended to excited states Yanai et al. (2005, 2015); Kottmann et al. (2015) and electric Sekino et al. (2008, 2012) and magnetic Jensen et al. (2016) linear response properties. In this approach all functions and operators, such as orbitals, densities and potential energy contributions to the Kohn–Sham Hamiltonian, are represented using Multiwavelets. Concerning potential energy terms, the external potential is obtained by projection ^{1}^{1}1For the singularity of the nuclear potential, a simple smoothing is employed, but its effect on the accuracy is easily controlled by a single parameter Harrison et al. (2004). onto the MW basis, the Hartree potential is computed through Poisson’s equation:
(6) 
and the XC potential is computed explicitly in the MW representation from the following expression:^{2}^{2}2Displayed as spinunpolarized for clarity. Extension to spinDFT is fairly straightforward.
(7) 
The partial derivatives of the XC kernel can be mapped pointwise though external XC libraries, Marques et al. (2012); Ekström et al. (2010) and the gradients are computed by the approach of Alpert et al. Alpert et al. (2002) For hybrid functionals, a fraction of the exact HartreeFock exchange contribution is included in the KS Hamiltonian (25% for PBE0). In our work we follow the method proposed by Yanai and coworkers, Yanai et al. (2004a) where the exchange operator is defined as:
(8) 
The above expression is again computed directly within the MW framework through repeated application of the Poisson operator to different orbital products.
While MWs are able to provide highaccuracy solution of integral equations in the form of Eq. (3), the same is not true for differential operators. In particular, highorder derivatives should be avoided in order to maintain accuracy in numerical algorithms. For this reason, we have found that the direct evaluation of the kinetic energy as a 2nd derivative of the wave function does not give the desired accuracy. Instead, we avoid the kinetic operator by computing the update to the eigenvalue directly:^{3}^{3}3This update is exact, provided that the orbital update comes directly from the application of the BSH operator defined by the previous (not necessarily exact) eigenvalue: . Generalizations can be made for multiple orbitals.
(9) 
where the ’s refer to differences between iterations and . In contrast, the gradients in the expression for the GGA potential in Eq. (7) have not been found to affect the accuracy, partly because of a slightly conservative over representation of the density grid,^{4}^{4}4The grid is constructed such that it holds both the density and its gradient within the requested accuracy. and partly because the XC energy is (by construction) only a small part of the total energy, thus reducing its relative accuracy requirement.
In this work, the MW calculations are performed with MRChem, mrc (2016) the GTO Feller (1996); Schuchardt et al. (2007) calculations with NWChem Valiev et al. (2010) and the NAO Delley (1990) calculations with FHIaims. Blum et al. (2009); Ren et al. (2012) APW+local orbital (APW+lo) calculations are performed with ELK. elk (2015) The exchangecorrelation functionals are calculated using the libxc Marques et al. (2012) library in case of NWChem and FHIaims, and the xcfun Ekström et al. (2010) library for MRChem.
The raw data of our study, as well as instructions for its reproducibility is available in the Supporting Information (SI).Jensen et al. (2017) Our test set comprises 211 molecules. In addition to the 147 systems from the G2/97 test set Schmider and Becke (1998) containing light elements up to the third row, it contains molecules with chemical elements that are underrepresented in the G2/97 test set (Be, Li, Mg, Al, F, Na, S and Cl) as well as 6 nonbonded systems. For most of the systems, the experimental structure obtained from the NIST Computational Chemistry Comparison and Benchmark Database NIS (2015) was employed. In the remaining cases, geometries have been optimized at the MP2 level of theory, using the largest Gaussian basis set (see SIJensen et al. (2017) for details).
We have considered four different basis sets each within NWChem and FHIaims, which include small ones intended for prerelaxations and energy differences between bonded structures ("light", 631+G**), production basis sets considered in most publications ("tight", “tier2” for FHIaims and augccpVDZ, augccpVTZ for Gaussian codes), as well the largest available basis set: “tier4” for FHIaims (“tier3” for H) and augccpV5Z for NWChem.^{5}^{5}5For Li, Be, Na and Mg the largest basis set is augccpVQZ and has therefore been employed. For the other elements, the corresponding 6Z basis is also available, but attempts to employ such a basis led to overcompleteness problems, often resulting in energies higher than the 5Z results.
An accurate, global resolutionofidentity approach ("RIV" in Ref.Ren et al. (2012)) is employed to evaluate the fourcenter Coulomb operator in hybridPBE0 in FHIaims. It is important to note that the NAO basis sets include a “minimal basis” of atomic radial functions determined for the same XC functional as used later in the threedimensional SCF calculations. This is standard practice in FHIaims for semilocal density functionals. For hybridPBE0, these radial functions are provided by linking FHIaims to the "atom_sphere" atomic solver code for spherically symmetric free atoms developed in the Goedecker group for several years.Goedecker and Maschke (1992)
The ground state energy of atoms from Hydrogen (=1) to Argon (=18) has been computed with the three chosen functionals. Our results are summarized in Figure 1. For all computational methods employed, the results of this section refer to the most accurate basis set employed (see previous section and SIJensen et al. (2017) for details). The top panel reports the LDASVWN5 values as absolute errors with respect to the reference values of the NIST Kotochigova et al. (1997) database for nonrelativistic, spin polarized, spherically symmetric atoms. As expected, MWs yield differences which are consistently below the requested accuracy of 1 Ha. The NAO and APW+lo approaches achieve average errors of 0.010.1 mHa and 0.11 mHa, respectively. GTOs are limited to around mHa accuracy. The GTO outliers (Li, Be, Na and Mg) have been computed with the augccpVQZ basis, because the augccpV5Z basis set is not available for these elements. Had 5Zquality functions been available for all elements, a more uniform error for GTO would have resulted along the series, but the overall picture would only improve slightly.
In the GGAPBE (middle) and hybridPBE0 (bottom) panels, all 18 atoms (both spherical and not), are included. The nonrelativistic, spinpolarized electronic density and the total energy of the ground state, computed using MRChem (converged within Ha) serves as the reference to which the other approaches are compared. For both functionals NWChem performs at the limit of chemical accuracy ( 1 mHa). The NAOs in FHIaims achieve 0.1 mHa or better, except for fluorine (0.3 mHa). For closedshell atoms, FHIaims is essentially exact because the exact radial functions of spherically symmetric, spinunpolarized atoms are included in the basis sets. For GTOs, we observe that the total energy error grows with the atomic number, . In contrast, the accuracy of NAOs is less affected by the nuclear charge, with errors generally below 0.1 mHa for the range examined here, irrespective of the choice of functional. For APW+lo, only the LDASVWN5 values are included in Figure 1: the corresponding GGAPBE and hybridPBE0 errors achieved in this work are above the threshold of 1e03 Ha (dashed line) and were not considered further because it is unclear how much they might be affected by implementationspecific aspects other than the basis set.
The total energies, atomization energies and dipole moments of the 211 molecules considered have been computed within the GGAPBE and hybridPBE0 functionals using MRChem with the highest affordable precision (below 1 Ha throughout). Fig. 2 reports the Mean Absolute Deviation (MAD), Root Mean Square Deviation (RMSD) and Maximum Absolute Deviation (maxAD) obtained for total energy (top panel), atomization energy (medium panel) and dipole moment (bottom panel) w.r.t MRChem for the GGAPBE and hybridPBE0 functionals, respectively.^{6}^{6}6Due to technical reasons in convergence, CHCHO was excluded from the PBE0 results, while CCH was excluded in both PBE and PBE0. For all the molecules, the correct groundstate spin multiplicity was specified.
Total energies are a measure of the accuracy achieved by each method/basis pair, whereas the atomization energies deserve special attention for their role in the development of density functionals, generally benchmarked against such thermodynamic values. However, as recently pointed out by Medvedev et al. Medvedev et al. (2017) the variational energy is not the optimal measure for the quality of the calculated electronic density, which influences numerous other observables. For this reason, we have included the dipole moment as a nonvariational quantity in our benchmarks (dipole errors are linear in the density error, whereas energies are quadratic). Dipole moments also serve as a verification that the different methods converge to the same electronic state and not to a nearby metastable configuration. Although the existence of multiple metastable SCF solutions in KohnSham DFT is well known, it is often not detected by users of electronic structure codes. The solution strategy, also employed in the present paper, is to probe different spin initializations of each molecule to identify the global minimum. In the present work, the correct identification has been validated by ensuring consistency of the dipole moment as well as the KS eigenvalue spectra produced by the three distinct electronic structure methods.
Several important conclusions can be drawn from the results obtained:

For total energies, FHIaims is able to reach more accurate results than NWChem with GTOs, for basis sets of comparable size (e.g. "tier4" vs. augccpV5Z).

For atomization energies, both NAOs and GTOs benefit from error cancellation to some extent. Such a cancellation is however much stronger for GTOs where the RMSD is lowered by a factor 48 in most cases, whereas for NAOs only by a factor of 1,52. In both cases the cancellation is more marked for the smallest bases. Despite the smaller cancellations, NAOs are still closer to the converged limit than GTOs, for comparable basis sets.

Due to the cumbersome convergence of periodic DFT codes with respect to the box size, we did not include APW+lo results for the entire test set of molecules. Nevertheless, for a small subset of molecules for which the limit of the box size was reached, we found atomization energies with errors of about 1 kcal/mol (see SIJensen et al. (2017)). Our experience suggests that it is technically challenging for APWbased codes to reach accuracies below 1 kcal/mol on atomization energies.
As a final remark, we stress that for a few atoms (Li, Be, Na, Mg), the augccpV5Z basis is not available, as previously mentioned in the atomic calculations part. Had it been available, GTOs might have yielded somewhat higher precision for the affected systems than in our benchmarks. However, considering the large size of our sample, the fact that only a few atoms in a molecule are affected, and the small improvement that can be inferred from the atomic calculations, our main conclusions still hold. On the other hand, such a de facto limitation of the availability of GTO basis sets illustrates how demanding it is to generate such basis sets. In contrast, MWs and NAOs are much less affected by such a limitation.
To the best of our knowledge, this work presents the most accurate atomization energies calculated to date, for a large benchmark set of molecules. We conclude that moderately sized GTO basis sets, frequently used in quantum chemistry applications, suffer from average total energy errors much larger than 10 kcal/mol, and while very large GTO basis sets yield the desired accuracy on average, there are still significant outliers. Furthermore, it may not always be feasible to employ such basis sets for systems much larger than those included in this study.
NAOs give much better accuracy even for moderately large bases (“tight” and beyond) since they can be constructed to possess the numerically correct behavior for a given XC functional, both in the nuclear as well as in the tail region. When feasible, APW+lobased calculations achieve errors around 1 mHa for total energies, and 1 kcal/mol for atomization energies. However, this level of convergence is difficult to reach for general molecular systems.
Our results show that the basis set error can dominate over errors arising from the choice of XC functional under many circumstances, in particular if some of the most advanced and accurate functionals are used. Our results set new standards in the verification and validation of electronic structure methods, and we expect them to be used to assess the accuracy of all future developments in Density Functional Theory methods.
Acknowledgements.
This research was partly supported by the Research Council of Norway through a Centre of Excellence Grant (Grant No. 179568/V30) with computing time provided by NOTUR (Grant No. NN4654K), and partly by the NCCR MARVEL, funded by the Swiss National Science Foundation (SNSF) with computing time provided by CSCS under project s707. S. Saha acknowledges support from the SNSF. J.A.F.L. acknowledges fruitful discussions with John Kay Dewhurst on the APW method.References
 Hohenberg and Kohn (1964) H. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
 Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
 Goerigk and Grimme (2010) L. Goerigk and S. Grimme, J. Chem. Theory Comput. 6, 107 (2010).
 Karton et al. (2008) A. Karton, A. Tarnopolsky, J.F. Lamere, G. C. Schatz, and J. M. L. Martin, J. Phys. Chem. A 112, 12868 (2008).
 Zhao and Truhlar (2005) Y. Zhao and D. G. Truhlar, J. Phys. Chem. A 109, 5656 (2005).
 Perdew et al. (1996) J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
 Perdew and Schmidt (2001) J. P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001).
 Li et al. (2015) C. Li, X. Zheng, A. J. Cohen, P. MoriSánchez, and W. Yang, Phys. Rev. Lett. 114, 053001 (2015).
 Sun et al. (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
 Su and Xu (2015) N. Q. Su and X. Xu, Int. J. Quantum Chem. 115, 589 (2015).
 Peverati and Truhlar (2014) R. Peverati and D. G. Truhlar, Philos. Trans. R. Soc. London A: Math., Phys. and Eng. Sci. 372, 20120476 (2014).
 Mardirossian and HeadGordon (2016) N. Mardirossian and M. HeadGordon, J. Chem. Theory Comput. 12, 4303 (2016), pMID: 27537680, http://dx.doi.org/10.1021/acs.jctc.6b00637 .
 Singh and Nordstrom (2006) D. J. Singh and L. Nordstrom, Planewaves, Pseudopotentials, and the LAPW method (Springer Science & Business Media, 2006).
 Goedecker and Maschke (1992) S. Goedecker and K. Maschke, Phy. Rev. A 45, 88 (1992).
 Willand et al. (2013) A. Willand, Y. O. Kvashnin, L. Genovese, Á. V. Mayagoitia, A. K. Deb, A. Sadeghi, T. Deutsch, and S. Goedecker, J. Chem. Phys. 138, 104109 (2013).
 Lejaeghere et al. (2016) K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, et al., Science 351, aad3000 (2016).
 Moncrieff and Wilson (2005) D. Moncrieff and S. Wilson, Int. J. Quantum Chem. 101, 363 (2005).
 Kotochigova et al. (1997) S. Kotochigova, Z. H. Levine, E. L. Shirley, M. D. Stiles, and C. W. Clark, Atomic reference data for electronic structure calculations (National Institute of Standards and Technology, Physics Laboratory, 1997).
 Ballone and Galli (1990) P. Ballone and G. Galli, Phys. Rev. B 42, 1112 (1990).
 Becke (1992) A. D. Becke, J. Chem. Phys. 96, 2155 (1992).
 Daykov et al. (2003) I. P. Daykov, T. A. Arias, and T. D. Engeness, Phys. Rev. Lett. 90 (2003), 10.1103/PhysRevLett.90.216402.
 NIS (2015) “NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Editor: Russell D. Johnson III,” See http://cccbdb.nist.gov/ (2015).
 Papas and Schaefer (2006) B. N. Papas and H. F. Schaefer, J. Mol. Struct.: THEOCHEM 768, 175 (2006).
 Harrison et al. (2003) R. J. Harrison, G. I. Fann, T. Yanai, and G. Beylkin, Comput. Sci.  ICCS 2003: Int. Conf. Proc. , 103 (2003).
 Harrison et al. (2004) R. J. Harrison, G. I. Fann, T. Yanai, Z. Gan, and G. Beylkin, J. Chem. Phys. 121, 11587 (2004).
 Yanai et al. (2004a) T. Yanai, G. I. Fann, Z. Gan, R. J. Harrison, and G. Beylkin, J. Chem. Phys. 121, 6680 (2004a).
 Yanai et al. (2004b) T. Yanai, G. I. Fann, Z. Gan, R. J. Harrison, and G. Beylkin, J. Chem. Phys. 121, 2866 (2004b).
 Vosko et al. (1980) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
 Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
 Paier et al. (2005) J. Paier, R. Hirschl, M. Marsman, and G. Kresse, J. Chem. Phys. 122, 234102 (2005).
 Medvedev et al. (2017) M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 355, 49 (2017).
 Pask and Sterne (2005) J. E. Pask and P. A. Sterne, Modelling and Simulation in Materials Science and Engineering 13, R71 (2005).
 Losilla and Sundholm (2012) S. A. Losilla and D. Sundholm, J. Chem. Phys. 136, 214104 (2012).
 Toivanen et al. (2015) E. A. Toivanen, S. A. Losilla, and D. Sundholm, Phys. Chem. Chem. Phys. 17, 31480 (2015).
 Parkkinen et al. (2017) P. Parkkinen, S. A. Losilla, E. Solala, E. A. Toivanen, W.H. Xu, and D. Sundholm, J. Chem. Theory and Comput. (2017), dOI:10.1021/acs.jctc.6b01207, http://dx.doi.org/10.1021/acs.jctc.6b01207 .
 Lin et al. (2012) L. Lin, J. Lu, L. Ying, and W. E, J. Comput. Phys. 231, 2140 (2012).
 Frediani and Sundholm (2015) L. Frediani and D. Sundholm, Phys. Chem. Chem. Phys. 17, 31357 (2015).
 Sundholm and Olsen (1994) D. Sundholm and J. Olsen, Phys. Rev. A 49, 3453 (1994).
 Alpert et al. (1993) B. K. Alpert, G. Beylkin, and R. Coifman, SIAM J. Sci. Stat. Comput. 14, 159 (1993).
 Alpert (1999) B. K. Alpert, SIAM J. Math. Analysis 24, 246 (1999).
 Beylkin and Mohlenkamp (2002) G. Beylkin and M. Mohlenkamp, Proc. Natl. Acad. Sci. USA 99, 10246 (2002).
 Beylkin et al. (2007) G. Beylkin, R. Cramer, G. Fann, and R. J. Harrison, Appl. Comput. Harmon. A. 23, 235 (2007).
 Fann et al. (2004) G. Fann, G. Beylkin, R. J. Harrison, and K. E. Jordan, IBM J. Research and Development 48, 161 (2004).
 Beylkin (2005) G. Beylkin, Appl. Comput. Harmon. A. 19, 17 (2005).
 Gines et al. (1998) D. Gines, G. Beylkin, and J. Dunn, Appl. Comput. Harmon. A. 5, 156 (1998).
 Beylkin et al. (2008) G. Beylkin, V. Cheruvu, and F. Perez, Appl. Comput. Harmon. A. 24, 354 (2008).
 Kalos (1962) M. H. Kalos, Phys. Rev. 128, 1791 (1962).
 Harrison (2004) R. J. Harrison, J. Comput. Chem. 25, 328 (2004).
 Yanai et al. (2005) T. Yanai, R. J. Harrison, and N. C. Handy, Mol. Phys. 103, 413 (2005).
 Yanai et al. (2015) T. Yanai, G. I. Fann, G. Beylkin, and R. J. Harrison, Phys. Chem. Chem. Phys. 17, 31405 (2015).
 Kottmann et al. (2015) J. S. Kottmann, S. Höfener, and F. A. Bischoff, Phys. Chem. Chem. Phys. 17, 31453 (2015).
 Sekino et al. (2008) H. Sekino, Y. Maeda, T. Yanai, and R. J. Harrison, J. Chem. Phys. 129, 034111 (2008).
 Sekino et al. (2012) H. Sekino, Y. Yokoi, and R. J. Harrison, J. Phys.: Conf. Series 352, 012014 (2012).
 Jensen et al. (2016) S. R. Jensen, T. Flå, D. Jonsson, R. S. Monstad, K. Ruud, and L. Frediani, Phys. Chem. Chem. Phys. 18, 21145 (2016).
 (55) For the singularity of the nuclear potential, a simple smoothing is employed, but its effect on the accuracy is easily controlled by a single parameter Harrison et al. (2004).
 (56) Displayed as spinunpolarized for clarity. Extension to spinDFT is fairly straightforward.
 Marques et al. (2012) M. A. Marques, M. J. Oliveira, and T. Burnus, Comput. Phys. Comm. 183, 2272 (2012).
 Ekström et al. (2010) U. Ekström, L. Visscher, R. Bast, A. J. Thorvaldsen, and K. Ruud, J. Chem. Theory Comput. 6, 1971 (2010).
 Alpert et al. (2002) B. K. Alpert, G. Beylkin, D. Gines, and L. Vozovoi, J. Comput. Phys. 182, 149 (2002).
 (60) This update is exact, provided that the orbital update comes directly from the application of the BSH operator defined by the previous (not necessarily exact) eigenvalue: . Generalizations can be made for multiple orbitals.
 (61) The grid is constructed such that it holds both the density and its gradient within the requested accuracy.
 mrc (2016) “Multiresolution chemistry (mrchem) program package.” See http://mrchemdoc.readthedocs.org/en/latest/ (2016).
 Feller (1996) D. Feller, J. Comput. Chem. 17, 1571 (1996).
 Schuchardt et al. (2007) K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, and T. L. Windus, J. Chem. Information and Modeling 47, 1045 (2007).
 Valiev et al. (2010) M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. van Dam, D. Wang, J. Nieplocha, E. Apra, T. Windus, and W. de Jong, Comput. Phys. Comm. 181, 1477 (2010).
 Delley (1990) B. Delley, J. Chem. Phys. 92, 508 (1990).
 Blum et al. (2009) V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Comput. Phys. Comm. 180, 2175 (2009).
 Ren et al. (2012) X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, and M. Scheffler, New Journal of Physics 14, 053020 (2012).
 elk (2015) “An allelectron fullpotential linearised augmentedplane wave (fplapw) code.” See http://http://elk.sourceforge.net/ (2015).
 Jensen et al. (2017) S. R. Jensen, S. Saha, J. A. FloresLivas, W. Huhn, V. Blum, S. Goedecker, and L. Frediani, “Ggapbe and hybridpbe0 energies and dipole moments with mrchem, fhiaims, nwchem and elk,” (2017).
 Schmider and Becke (1998) H. L. Schmider and A. D. Becke, J. Chem. Phys. 108, 9624 (1998).
 (72) For Li, Be, Na and Mg the largest basis set is augccpVQZ and has therefore been employed. For the other elements, the corresponding 6Z basis is also available, but attempts to employ such a basis led to overcompleteness problems, often resulting in energies higher than the 5Z results.
 (73) Due to technical reasons in convergence, CHCHO was excluded from the PBE0 results, while CCH was excluded in both PBE and PBE0.
 Bak et al. (2000) K. L. Bak, J. Gauss, T. Helgaker, P. Jørgensen, and J. Olsen, Chem. Phys. Lett. 319, 563 (2000).
 Ha
 Hartree
 CBS
 complete basis set
 MRA
 Multiresolution Analysis
 MW
 Multiwavelet
 SI
 Supporting Information
 GTO
 Gaussiantype orbital
 BSH
 BoundState Helmholtz
 SCF
 SelfConsistent Field
 HF
 Hartree–Fock
 KS
 Kohn–Sham
 XC
 ExchangeCorrelation
 DFT
 Density Functional Theory
 GGA
 Generalized Gradient Approximation
 LDA
 Local Density Approximation
 KAIN
 KrylovAccelerated Inexact Newton
 AO
 Atomic Orbital
 MO
 Molecular Orbital
 LCAO
 Linear Combination of Atomic Orbitals
 NAO
 numeric atomcentered orbital
 LAPW
 linearised augmented plane wave
 APW
 augmented plane wave
 APW+lo
 APW+local orbital
 MAD
 Mean Absolute Deviation
 MD
 Mean Deviation
 maxAD
 Maximum Absolute Deviation
 minAD
 Minimum Absolute Deviation
 RMSD
 Root Mean Square Deviation
 SD
 Standard Deviation
 MBPT
 Manybody perturbation theory
 MP2
 2ndorder MøllerPlesset
 CC
 Coupled Cluster
 CI
 Configuration Interaction
 CCSD
 CoupledCluster Singles Doubles
 MCSCF
 MultiConfigurational SelfConsistent Field