The Elephant in the Room of Density Functional Theory Calculations

# The Elephant in the Room of Density Functional Theory Calculations

## 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 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 10 kcal/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.

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-GGAsSun 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 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 (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 milli-Hartree (mHa), and even 10 decimal places for micro-Hartree (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 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 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 , such as the Legendre or the Interpolating polynomials, re-scaled 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:

 V0k⊂V1k⊂⋯⊂Vnk⊂⋯⊂L2. (1)

The wavelet spaces are simply the orthogonal complement of two subsequent scaling spaces and :

 Wnk⊕Vnk=Vn+1k,Wnk⊥Vnk. (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 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:

 g(r)=^Gμf(r)=∫Gμ(r−r′)f(r′)dr′. (3)

Among such operators are the Poisson () and the Bound-State Helmholtz (BSH) () kernels:

 Gμ(r−r′)=e−μ|r−r′|4π|r−r|. (4)

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)

 φi=−2^Gμi^Vφi, (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 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 1 onto the MW basis, the Hartree potential is computed through Poisson’s equation:

 VHartree(r)=∫ρ(r′)4π|r−r′|dr′, (6)

and the XC potential is computed explicitly in the MW representation from the following expression:2

 VXC(r)=∂fXC∂ρ−2∇∂fXC∂|∇ρ|2⋅∇ρ. (7)

The partial derivatives of the XC kernel 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:

 ^Kf(r)=∑jφj(r)∫φj(r′)f(r′)4π|r−r′|dr′. (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 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:3

 Δϵn=⟨φn+1|^Vn|Δφn⟩⟨φn+1|φn+1⟩+⟨φn+1|Δ^Vn|φn+1⟩⟨φn+1|φn+1⟩, (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 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 MRChemmrc (2016) the GTO Feller (1996); Schuchardt et al. (2007) calculations with NWChem Valiev et al. (2010) and the NAO Delley (1990) calculations with FHI-aimsBlum et al. (2009); Ren et al. (2012) APW+local orbital (APW+lo) calculations are performed with ELKelk (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.5

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 (=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 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 Ha. The NAO and APW+lo approaches achieve average errors of 0.01-0.1 mHa and 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 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 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, . 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 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 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.6 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:

1. 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. 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. 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. 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. 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.

###### 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.
Ha
Hartree
CBS
complete basis set
MRA
Multiresolution Analysis
MW
Multiwavelet
SI
Supporting Information
GTO
Gaussian-type orbital
BSH
Bound-State Helmholtz
SCF
Self-Consistent Field
HF
Hartree–Fock
KS
Kohn–Sham
XC
Exchange-Correlation
DFT
Density Functional Theory
GGA
LDA
Local Density Approximation
KAIN
Krylov-Accelerated Inexact Newton
AO
Atomic Orbital
MO
Molecular Orbital
LCAO
Linear Combination of Atomic Orbitals
NAO
numeric atom-centered orbital
LAPW
linearised augmented plane wave
APW
augmented plane wave
APW+lo
APW+local orbital
Mean Absolute Deviation
MD
Mean Deviation
Maximum Absolute Deviation
Minimum Absolute Deviation
RMSD
Root Mean Square Deviation
SD
Standard Deviation
MBPT
Many-body perturbation theory
MP2
2nd-order Møller-Plesset
CC
Coupled Cluster
CI
Configuration Interaction
CCSD
Coupled-Cluster Singles Doubles
MCSCF
Multi-Configurational Self-Consistent Field

### Footnotes

1. 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).
2. Displayed as spin-unpolarized for clarity. Extension to spin-DFT is fairly straightforward.
3. 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.
4. The grid is constructed such that it holds both the density and its gradient within the requested accuracy.
5. For 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.
6. Due to technical reasons in convergence, CHCHO was excluded from the PBE0 results, while CCH was excluded in both PBE and PBE0.

### References

1. H. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
2. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
3. L. Goerigk and S. Grimme, J. Chem. Theory Comput. 6, 107 (2010).
4. A. Karton, A. Tarnopolsky, J.-F. Lamere, G. C. Schatz,  and J. M. L. Martin, J. Phys. Chem. A 112, 12868 (2008).
5. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A 109, 5656 (2005).
6. J. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
7. J. P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001).
8. C. Li, X. Zheng, A. J. Cohen, P. Mori-Sánchez,  and W. Yang, Phys. Rev. Lett. 114, 053001 (2015).
9. J. Sun, A. Ruzsinszky,  and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
10. N. Q. Su and X. Xu, Int. J. Quantum Chem. 115, 589 (2015).
11. R. Peverati and D. G. Truhlar, Philos. Trans. R. Soc. London A: Math., Phys. and Eng. Sci. 372, 20120476 (2014).
12. N. Mardirossian and M. Head-Gordon, J. Chem. Theory Comput. 12, 4303 (2016), pMID: 27537680, http://dx.doi.org/10.1021/acs.jctc.6b00637 .
13. D. J. Singh and L. Nordstrom, Planewaves, Pseudopotentials, and the LAPW method (Springer Science & Business Media, 2006).
14. S. Goedecker and K. Maschke, Phy. Rev. A 45, 88 (1992).
15. 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).
16. 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).
17. D. Moncrieff and S. Wilson, Int. J. Quantum Chem. 101, 363 (2005).
18. 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).
19. P. Ballone and G. Galli, Phys. Rev. B 42, 1112 (1990).
20. A. D. Becke, J. Chem. Phys. 96, 2155 (1992).
21. I. P. Daykov, T. A. Arias,  and T. D. Engeness, Phys. Rev. Lett. 90 (2003), 10.1103/PhysRevLett.90.216402.
22. “NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Editor: Russell D. Johnson III,” See http://cccbdb.nist.gov/ (2015).
23. B. N. Papas and H. F. Schaefer, J. Mol. Struct.: THEOCHEM 768, 175 (2006).
24. R. J. Harrison, G. I. Fann, T. Yanai,  and G. Beylkin, Comput. Sci. - ICCS 2003: Int. Conf. Proc. , 103 (2003).
25. R. J. Harrison, G. I. Fann, T. Yanai, Z. Gan,  and G. Beylkin, J. Chem. Phys. 121, 11587 (2004).
26. T. Yanai, G. I. Fann, Z. Gan, R. J. Harrison,  and G. Beylkin, J. Chem. Phys. 121, 6680 (2004a).
27. T. Yanai, G. I. Fann, Z. Gan, R. J. Harrison,  and G. Beylkin, J. Chem. Phys. 121, 2866 (2004b).
28. S. H. Vosko, L. Wilk,  and M. Nusair, Can. J. Phys. 58, 1200 (1980).
29. C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
30. J. Paier, R. Hirschl, M. Marsman,  and G. Kresse, J. Chem. Phys. 122, 234102 (2005).
31. M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew,  and K. A. Lyssenko, Science 355, 49 (2017).
32. J. E. Pask and P. A. Sterne, Modelling and Simulation in Materials Science and Engineering 13, R71 (2005).
33. S. A. Losilla and D. Sundholm, J. Chem. Phys. 136, 214104 (2012).
34. E. A. Toivanen, S. A. Losilla,  and D. Sundholm, Phys. Chem. Chem. Phys. 17, 31480 (2015).
35. 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 .
36. L. Lin, J. Lu, L. Ying,  and W. E, J. Comput. Phys. 231, 2140 (2012).
37. L. Frediani and D. Sundholm, Phys. Chem. Chem. Phys. 17, 31357 (2015).
38. D. Sundholm and J. Olsen, Phys. Rev. A 49, 3453 (1994).
39. B. K. Alpert, G. Beylkin,  and R. Coifman, SIAM J. Sci. Stat. Comput. 14, 159 (1993).
40. B. K. Alpert, SIAM J. Math. Analysis 24, 246 (1999).
41. G. Beylkin and M. Mohlenkamp, Proc. Natl. Acad. Sci. USA 99, 10246 (2002).
42. G. Beylkin, R. Cramer, G. Fann,  and R. J. Harrison, Appl. Comput. Harmon. A. 23, 235 (2007).
43. G. Fann, G. Beylkin, R. J. Harrison,  and K. E. Jordan, IBM J. Research and Development 48, 161 (2004).
44. G. Beylkin, Appl. Comput. Harmon. A. 19, 17 (2005).
45. D. Gines, G. Beylkin,  and J. Dunn, Appl. Comput. Harmon. A. 5, 156 (1998).
46. G. Beylkin, V. Cheruvu,  and F. Perez, Appl. Comput. Harmon. A. 24, 354 (2008).
47. M. H. Kalos, Phys. Rev. 128, 1791 (1962).
48. R. J. Harrison, J. Comput. Chem. 25, 328 (2004).
49. T. Yanai, R. J. Harrison,  and N. C. Handy, Mol. Phys. 103, 413 (2005).
50. T. Yanai, G. I. Fann, G. Beylkin,  and R. J. Harrison, Phys. Chem. Chem. Phys. 17, 31405 (2015).
51. J. S. Kottmann, S. Höfener,  and F. A. Bischoff, Phys. Chem. Chem. Phys. 17, 31453 (2015).
52. H. Sekino, Y. Maeda, T. Yanai,  and R. J. Harrison, J. Chem. Phys. 129, 034111 (2008).
53. H. Sekino, Y. Yokoi,  and R. J. Harrison, J. Phys.: Conf. Series 352, 012014 (2012).
54. 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 spin-unpolarized for clarity. Extension to spin-DFT is fairly straightforward.
57. M. A. Marques, M. J. Oliveira,  and T. Burnus, Comput. Phys. Comm. 183, 2272 (2012).
58. U. Ekström, L. Visscher, R. Bast, A. J. Thorvaldsen,  and K. Ruud, J. Chem. Theory Comput. 6, 1971 (2010).
59. 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.
62. “Multiresolution chemistry (mrchem) program package.” See http://mrchemdoc.readthedocs.org/en/latest/ (2016).
63. D. Feller, J. Comput. Chem. 17, 1571 (1996).
64. 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).
65. 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).
66. B. Delley, J. Chem. Phys. 92, 508 (1990).
67. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter,  and M. Scheffler, Comput. Phys. Comm. 180, 2175 (2009).
68. X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter,  and M. Scheffler, New Journal of Physics 14, 053020 (2012).
69. “An all-electron full-potential linearised augmented-plane wave (fp-lapw) code.” See http://http://elk.sourceforge.net/ (2015).
70. S. R. Jensen, S. Saha, J. A. Flores-Livas, W. Huhn, V. Blum, S. Goedecker,  and L. Frediani, “Gga-pbe and hybrid-pbe0 energies and dipole moments with mrchem, fhi-aims, nwchem and elk,”  (2017).
71. 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 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.
73. Due to technical reasons in convergence, CHCHO was excluded from the PBE0 results, while CCH was excluded in both PBE and PBE0.
74. K. L. Bak, J. Gauss, T. Helgaker, P. Jørgensen,  and J. Olsen, Chem. Phys. Lett. 319, 563 (2000).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters