The Elephant in the Room of Density Functional Theory Calculations

TL;DR

This study uses multiwavelet methods to generate highly accurate reference energies for molecules, enabling a detailed assessment of the accuracy of various basis sets in density functional theory calculations.

Contribution

It provides the first systematic quantification of basis set errors in DFT using quasi-exact reference energies, highlighting the accuracy limits of common basis sets.

Findings

NAOs achieve chemical accuracy with typical basis sizes

GTOs require larger basis sets for similar accuracy

Systematic error cancellations improve atomization energy accuracy

Abstract

Using multiwavelets, we have obtained total energies and corresponding atomization energies for the GGA-PBE and hybrid-PBE0 density functionals for a test set of 211 molecules with an unprecedented and guaranteed $\mu$Hartree accuracy. These quasi-exact references allow us to quantify the accuracy of standard all-electron basis sets that are believed to be highly accurate for molecules, such as Gaussian-type orbitals (GTOs), all-electron numeric atom-centered orbitals (NAOs) and full-potential augmented plane wave (APW) methods. We show that NAOs are able to achieve the so-called chemical accuracy (1 kcal/mol) for the typical basis set sizes used in applications, for both total and atomization energies. For GTOs, a triple-zeta quality basis has mean errors of ~10kcal/mol in total energies, while chemical accuracy is almost reached for a quintuple-zeta 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.

Equations

$$V_k^0\subset V_k^1\subset \cdots \subset V_k^n\subset \cdots \subset L^2.$$

$$W_k^n\oplus V_k^n=V_k^{n+1},W_k^n\perp V_k^n.$$

$$g(r)={\hat{G}}_{\mu }f(r)=\int G_{\mu }(r-r^{\prime })f(r^{\prime })d{r}^{\prime }.$$

$$G_{\mu }(r-r^{\prime })=\frac{e^{-\mu |r-r^{\prime }|}}{4\pi |r-r|}.$$

$${\phi }_i=-2{\hat{G}}_{{\mu }_i}\hat{V}{\phi }_i,$$

$$V_{\text{Hartree}}(r)=\int \frac{\rho (r^{\prime })}{4\pi |r-r^{\prime }|}d{r}^{\prime },$$

$$V_{\text{XC}}(r)=\frac{\partial f_{\text{XC}}}{\partial \rho }-2\nabla \frac{\partial f_{\text{XC}}}{\partial |\nabla \rho {|}^2}\cdot \nabla \rho .$$

$$\hat{K}f(r)=\sum _j{\phi }_j(r)\int \frac{{\phi }_j(r^{\prime })f(r^{\prime })}{4\pi |r-r^{\prime }|}\mathrm{d}{r}^{\prime }.$$

$$\Delta {ϵ}^n=\frac{⟨{\phi }^{n+1}|{\hat{V}}^n|\Delta {\phi }^n⟩}{⟨{\phi }^{n+1}|{\phi }^{n+1}⟩}+\frac{⟨{\phi }^{n+1}|\Delta {\hat{V}}^n|{\phi }^{n+1}⟩}{⟨{\phi }^{n+1}|{\phi }^{n+1}⟩},$$

Full text

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 Coupled-Cluster 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 Exchange-Correlation (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) meta-GGAs, 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 main-group elements as well as transition metals and solids. Peverati and Truhlar (2014); Mardirossian and Head-Gordon (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 GGA-PBEPerdew et al. (1996) calculated equations of state for 71 elemental crystals from 15 different widely used DFT codes, employing both all-electron 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 basis-set choice: all-electron numeric atom-centered 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 Gaussian-type orbital (GTO) and NAO basis sets: both bases cannot be systematically enlarged to achieve completeness in the $L^{2}$ 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 self-consistent three-dimensional 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 (SiCl4) in our data set is less than 1 Ha, but it is computed as a difference of energies in the order of $\sim$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 $\sim$ milli-Hartree (mHa), and even 10 decimal places for micro-Hartree ($\mu$Ha) accuracy.

For isolated atoms, and using the appropriate numerical techniques, the associated many-particle 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 transform-based 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 non-orthogonality, which leads to inevitable algebraic ill-conditioning 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 $Z/|r-R|$ potential is replaced by a smoother potential that retains approximately the same physical properties of the all-electron 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 all-electron calculations. With this approach, it is now possible to achieve all-electron energies with arbitrarily small errors.

In the present work, we use MWs to obtain error bars of less than a $\mu$ 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, LDA-SVWN5, Vosko et al. (1980) GGA-PBE, Perdew et al. (1996) as well as hybrid-PBE0. 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 quasi-exact 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)

Real-space 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 so-called 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 multi-scale 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 non-uniform 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 $k$, such as the Legendre or the Interpolating polynomials, re-scaled and orthonormalized on the unit interval $[0,1]$. Then, an orthonormal scaling basis at refinement level $2^{-n}$ is constructed by dilation and translation of the original basis functions $\phi^{n}_{i,l}(x)=2^{n/2}\phi_{i}(2^{n}x-l)$, where $\phi^{n}_{i,l}$ is the $i$-th polynomial in the interval $[l/2^{n},(l+1)/2^{n})]$ at scale $n$. The set of scaling functions on all $2^{n}$ translations at scale $n$ defines the scaling space $V_{k}^{n}$, and in this way a ladder of spaces is constructed such that the complete $L^{2}$ limit can be approached in a systematic manner:

[TABLE]

The wavelet spaces $W^{n}_{k}$ are simply the orthogonal complement of two subsequent scaling spaces $V_{k}^{n}$ and $V_{k}^{n+1}$:

[TABLE]

Completeness in the $L^{2}$ sense can be achieved both by increasing the polynomial order (larger $k$) of the basis and by increasing the refinement in the ladder of spaces (larger $n$).

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 narrow-banded 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 tensor-product 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 multi-dimensional operator as a product of one-dimensional 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 non-standard 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, non-standard form of operators) made the accurate application of several important convolution operators efficient in three dimensions:

[TABLE]

Among such operators are the Poisson ($\mu=0$) and the Bound-State Helmholtz (BSH) ($\mu^{2}>0$) kernels:

[TABLE]

This mathematical framework was introduced to the computational chemistry community in the mid 2000’s by Harrison and co-workers. 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)

[TABLE]

where the $\varphi_{i}$’s are the KS orbitals and the potential operator $\hat{V}$ includes external (nuclear), Hartree, exchange and correlation contributions, while the kinetic operator and the orbital energies are included in the BSH operator as 2\hat{G}_{\mu_{i}}=\big{(}\hat{T}-\epsilon_{i}\big{)}^{-1}, with $\mu_{i}^{2}=-2\epsilon_{i}$. The ordinary KS equations $\hat{H}\varphi_{i}=\epsilon_{i}\varphi_{i}$, where $\hat{H}=\hat{T}+\hat{V}$ is the KS Hamiltonian, can be recovered by recalling that $(\nabla^{2}-\mu^{2})G_{\mu}(r-r^{\prime})=-\delta(r-r^{\prime})$. Such an integral formalism, when combined with a Krylov subspace accelerator, Harrison (2004) leads to fast and robust convergence of the fix-point iteration of Eq. (5).

The integral formulation, in combination with MWs, provides unprecedented accuracy for all-electron 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 111For 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:

[TABLE]

and the XC potential is computed explicitly in the MW representation from the following expression:222Displayed as spin-unpolarized for clarity. Extension to spin-DFT is fairly straightforward.

[TABLE]

The partial derivatives of the XC kernel $f_{XC}$ can be mapped point-wise 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 Hartree-Fock 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:

[TABLE]

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 high-accuracy solution of integral equations in the form of Eq. (3), the same is not true for differential operators. In particular, high-order 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:333This 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: $\Delta\varphi^{n}=-2\hat{G}^{n}_{\mu}\left[\hat{V}^{n}\varphi^{n}\right]-\varphi^{n}$. Generalizations can be made for multiple orbitals.

[TABLE]

where the $\Delta$’s refer to differences between iterations $n$ and $n+1$. 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,444The 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 FHI-aims. Blum et al. (2009); Ren et al. (2012) APW+local orbital (APW+lo) calculations are performed with ELK. elk (2015) The exchange-correlation functionals are calculated using the libxc Marques et al. (2012) library in case of NWChem and FHI-aims, 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 non-bonded 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 FHI-aims, which include small ones intended for prerelaxations and energy differences between bonded structures ("light", 6-31+G**), production basis sets considered in most publications ("tight", “tier2” for FHI-aims and aug-cc-pVDZ, aug-cc-pVTZ for Gaussian codes), as well the largest available basis set: “tier4” for FHI-aims (“tier3” for H) and aug-cc-pV5Z for NWChem.555For Li, Be, Na and Mg the largest basis set is aug-cc-pVQZ 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 resolution-of-identity approach ("RI-V" in Ref.Ren et al. (2012)) is employed to evaluate the four-center Coulomb operator in hybrid-PBE0 in FHI-aims. 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 three-dimensional SCF calculations. This is standard practice in FHI-aims for semilocal density functionals. For hybrid-PBE0, these radial functions are provided by linking FHI-aims 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 ($Z$=1) to Argon ($Z$=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 LDA-SVWN5 values as absolute errors with respect to the reference values of the NIST Kotochigova et al. (1997) database for non-relativistic, spin polarized, spherically symmetric atoms. As expected, MWs yield differences which are consistently below the requested accuracy of 1 $\mu$Ha. The NAO and APW+lo approaches achieve average errors of $\sim$0.01-0.1 mHa and $\sim$0.1-1 mHa, respectively. GTOs are limited to around mHa accuracy. The GTO outliers (Li, Be, Na and Mg) have been computed with the aug-cc-pVQZ basis, because the aug-cc-pV5Z basis set is not available for these elements. Had 5Z-quality 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 GGA-PBE (middle) and hybrid-PBE0 (bottom) panels, all 18 atoms (both spherical and not), are included. The non-relativistic, spin-polarized electronic density and the total energy of the ground state, computed using MRChem (converged within $\mu$Ha) serves as the reference to which the other approaches are compared. For both functionals NWChem performs at the limit of chemical accuracy ($\sim$ 1 mHa). The NAOs in FHI-aims achieve 0.1 mHa or better, except for fluorine (0.3 mHa). For closed-shell atoms, FHI-aims is essentially exact because the exact radial functions of spherically symmetric, spin-unpolarized atoms are included in the basis sets. For GTOs, we observe that the total energy error grows with the atomic number, $Z$. In contrast, the accuracy of NAOs is less affected by the nuclear charge, with errors generally below 0.1 mHa for the $Z$ range examined here, irrespective of the choice of functional. For APW+lo, only the LDA-SVWN5 values are included in Figure 1: the corresponding GGA-PBE and hybrid-PBE0 errors achieved in this work are above the threshold of 1e-03 Ha (dashed line) and were not considered further because it is unclear how much they might be affected by implementation-specific aspects other than the basis set.

The total energies, atomization energies and dipole moments of the 211 molecules considered have been computed within the GGA-PBE and hybrid-PBE0 functionals using MRChem with the highest affordable precision (below 1 $\mu$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 GGA-PBE and hybrid-PBE0 functionals, respectively.666Due to technical reasons in convergence, CH3CH2O was excluded from the PBE0 results, while CCH was excluded in both PBE and PBE0. For all the molecules, the correct ground-state 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 non-variational 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 Kohn-Sham 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, FHI-aims is able to reach more accurate results than NWChem with GTOs, for basis sets of comparable size (e.g. "tier4" vs. aug-cc-pV5Z). 2. 2.

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 4-8 in most cases, whereas for NAOs only by a factor of 1,5-2. 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. 3. 3.

The two functionals considered (GGA-PBE and hybrid-PBE0) yield very similar results, and we therefore assume that our conclusions concerning the accuracy of the different approaches (NAOs and GTOs) will hold also for other functionals of the same type. 4. 4.

Dipole moments can be considered accurate if deviations are below 0.01 Debye. Bak et al. (2000); NIS (2015) Only the largest basis sets in NAO and GTO used in our calculations achieve this target on average, but even such basis sets have outliers with errors close to 0.1 Debye. 5. 5.

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 APW-based 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 aug-cc-pV5Z 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+lo-based 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.