Meta-GGA Density Functional Calculations on Atoms with Spherically Symmetric Densities in the Finite Element Formalism

Density functional calculations on atoms are often used for determining accurate initial guesses as well as generating various types of pseudopotential approximations and efficient atomic-orbital basis sets for polyatomic calculations. To reach the best accuracy for these purposes, the atomic calculations should employ the same density functional as the polyatomic calculation. Atomic density functional calculations are typically carried out employing spherically symmetric densities, corresponding to the use of fractional orbital occupations. We have described their implementation for density functional approximations (DFAs) belonging to the local density approximation (LDA) and generalized gradient approximation (GGA) levels of theory as well as Hartree–Fock (HF) and range-separated exact exchange [Lehtola, S. Phys. Rev. A2020, 101, 012516]. In this work, we describe the extension to meta-GGA functionals using the generalized Kohn–Sham scheme, in which the energy is minimized with respect to the orbitals, which in turn are expanded in the finite element formalism with high-order numerical basis functions. Furnished with the new implementation, we continue our recent work on the numerical well-behavedness of recent meta-GGA functionals [Lehtola, S.; Marques, M. A. L. J. Chem. Phys.2022, 157, 174114]. We pursue complete basis set (CBS) limit energies for recent density functionals and find many to be ill-behaved for the Li and Na atoms. We report basis set truncation errors (BSTEs) of some commonly used Gaussian basis sets for these density functionals and find the BSTEs to be strongly functional dependent. We also discuss the importance of density thresholding in DFAs and find that all of the functionals studied in this work yield total energies converged to 0.1 μEh when densities smaller than 10–11a0–3 are screened out.


INTRODUCTION
Atoms are interesting for fundamental quantum chemistry, as they form the simplest bound many-electron systems. The key aspect of the electronic structure of atoms is their shell structure, which arises from the significant amount of symmetry inherent in the Coulomb problem. Importantly for chemistry, the shell structure of atoms is preserved to a large extent also in polyatomic systems, because the nuclear Coulomb potential dominates close to the nucleus, V(r) = −Z/r → − ∞ when r → 0, and thus the innermost electronic orbitals turn out to be insensitive to changes in the chemical environment. This feature arguably makes atomic calculations the keystone of electronic structure calculations: the near-constant nature of the shell structure is the assumption made in most computational approaches in the electronic structure theory of polyatomic systems, as we will discuss in the following.
Quantum chemical calculations on polyatomic systems invariably start from the solution of a self-consistent field (SCF) problem. 1 The iterative solution of the SCF problem requires an initial guess for the electron density, or the electronic orbitals. The best initial guesses are those that correctly reproduce the shell structure of atoms; 2 good alternatives include the superposition of atomic densities 3,4 (SAD) guess, the superposition of atomic potentials 2,5 (SAP) guess, as well as a parameter-free extended Huckel guess 6 that similarly can also be derived from atomic calculations. 2,7 Also the available numerical approaches used to carry out the polyatomic calculation are tightly connected to atomic electronic structure. The dominant basis set for electronic structure calculations in the solid state is plane waves; however, such calculations invariably employ pseudopotentials 8 or the projector augmented wave (PAW) method 9 that eliminate the need for an explicit description for the chemically inactive core electrons. Pseudopotentials and PAW setups are again derived from atomic calculations. All-electron plane-wave calculations have only extremely recently been shown to be feasible through the use of a regularized nuclear Coulomb potential, 10,11 but such calculations will likely be reserved for benchmark purposes due to their high computational cost.
In contrast, molecular quantum chemical calculations almost invariably employ the linear combination of atomic orbitals (LCAO) approach, in which the ith molecular orbital of spin σ, ψ iσ (r), is expressed as an expansion (2) which are composed of radial functions R nl labeled by the primary quantum number n and angular quantum number l, combined with spherical harmonics Y lm from m = −l to m = l in either the real or complex form, C αi σ being the corresponding expansion coefficients. Two major advantages of LCAO calculations are that (i) all electrons can be explicitly modeled, facilitating access to electronic core properties, for instance, and that (ii) the basis set truncation errors turn out to be systematic in many cases and cancel out in chemically relevant energy differences, affording near-quantitative accuracy with basis sets of modest size. 12 Despite their name, the (radial) AOs used in the LCAO expansion do not have to represent actual atomic orbitals, that is, physical one-particle states of an actual atoms. Instead, several types of radial functions such as Gaussian-type orbitals (GTOs) and Slater-type orbitals (STOs) may be and are commonly used in eq 2; see ref 12 for discussion. However, actual atomic orbitals are in principle the best option. These AOs can be solved by a fully numerical approach, 12 yielding so-called numerical atomic orbitals (NAOs). NAOs are an especially powerful basis set for electronic structure calculations: the minimal NAO basis is already exact for noninteracting atoms in SCF calculations. 12,13 Thanks to this exactness, the issues with basis set superposition error 14 (BSSE) that complicate the determination of reliable molecular geometries in calculations with GTOs or STOs are less of an issue in calculations employing NAOs, as "borrowing" basis functions from other atoms does not lead to improvement of atomic energies. Combined with carefully formed basis sets, 15 NAO methods have been shown to afford an excellent level of accuracy compared to fully numerical calculations. 16 A number of programs relying on NAOs in either all-electron or pseudopotential form have been published and are in active use by the community. 15,17−22 Because the chemical bonding situation of an atom�and the related deformation of the atom's electron density�is not known a priori in a molecule, atomic starting densities, starting potentials, pseudopotentials and NAOs are typically computed using spherically symmetric densities, achieved by fractional occupations of the atomic orbitals, as this ensures that all bonding situations are described equally, on the same footing. While such an approximation may sound coarse compared to the behavior of the real atom, the approximation does yield a correct shell structure, and thus offers a simple and sensible starting point for more sophisticated calculations, such as embedding the atom in a polyatomic calculation.
Still, to reach their best accuracy, atomic starting densities, starting potentials, pseudopotentials and NAOs should be determined for the same exchange-correlation functional that is employed in the polyatomic calculation. Because different functionals lead to different orbitals, the use of inconsistent NAOs may lead to the resurrection of issues with BSSE, for instance, as the sought-for exactness property is not satisfied in such a case.
However, meta-GGA functionals are not yet supported by many atomic solvers, especially when exact exchange is also included in the functional. 12 Although some programs have already been extended for meta-GGA functionals, 23−25 the approaches are often not fully self-consistent. For example, the work of Sun et al. 23 appears to have used the same PAWs for all functionals, while the recent work of Doumont et al. 26 is likewise unable to describe core electrons and generate atomic basis functions self-consistently with meta-GGA functionals, being limited to the GGA level for full self-consistency. Many other programs still lack support for fully self-consistent meta-GGAs. We will show in this work that the self-consistent implementation of meta-GGAs for atoms with spherical symmetry is not much more demanding than that of GGA functionals.
A major motivation of this work are repeated queries for reliable atomic reference data from developers of other fully numerical approaches. 27−29 Implementing any new computational method or algorithm requires being able to test whether the new implementation is correct, and the verification of any new atomic implementation for initial guesses, pseudopotentials or NAOs thereby requires access to reliable, high-quality reference data. Although the National Institute of Standards and Technology (NIST) hosts an atomic structure database, 30 its content is limited to calculations performed with the local density approximation (LDA). 31,32 While thorough data sets on some LDA and generalized gradient approximation (GGA) functionals can be found in the literature, 33,34 sub-μE h accurate reference energies for atoms with fractional occupations and meta-GGA functionals have not been published to this author's best knowledge. A key goal of this work is to provide such highly reliable reference values for use in verifying other implementations.
A further motivation of this work is the need to characterize and study the numerical behavior of meta-GGA functionals. The present author is a long-time developer of the Libxc library of density functionals, 35 which is used at present by some 40 electronic structure programs. We have recently thoroughly examined the numerical behavior of all the density functionals in Libxc at fixed atomic electron densities. 36 However, as discussed in ref 36, the ultimate test for the numerical stability of density functionals is the determination of complete basis set (CBS) limit energies in fully numerical calculations, as this requires accurate evaluation of the total energy in a sequence of numerical basis sets of increasing size. Fully numerical calculations are a demanding test of the well-behavedness of density functional approximations (DFAs). The ability to run fully numerical calculations with various types of DFAs is a great boon for the development of novel DFAs, as well as for their reliable implementation in Libxc.
For all of the above reasons, it would be appealing to be able to run self-consistent calculations with meta-GGA functionals quickly and reliably in extended basis sets. The finite element method (FEM) offers an attractive solution for determining reliable NAOs and total energies. FEM affords a variational approach to the CBS limit in fully numerical calculations, 12,37 and atomic Hartree−Fock ground-state energies converge extremely rapidly to the CBS limit with respect to the size of the radial basis set, when high-order numerical basis sets are used. 37, 38 We have previously published a general atomic solver in HelFEM 34,37,44 that is able to handle meta-GGA functionals including global hybrids with modern finite element approaches and verified it against established Gaussian-basis approaches in ref 37. We have also recently examined the role of the finite element shape functions in meta-GGA calculations, and found that the requirements for the numerical basis set are similar for HF as well as for density functional calculations with LDA, GGA and meta-GGA functionals, and that Lagrange and Hermite interpolating polynomials can either be used to pursue the CBS limit for τ-dependent meta-GGA functionals. 38 However, the solver described in ref 37 targets general wave functions, which may exhibit symmetry breaking, while most NAO approaches assume spherically symmetric orbitals with fractional occupations. We have recently discussed the extension of the FEM approach to the case of spherical symmetry with LDA and GGA orbitals as well as range-separated hybrids in ref 34, and reported nonrelativistic Hartree−Fock ground states for spin-unrestricted and spin-restricted calculations for H−Og (1 ≤ Z ≤ 118) in refs 5 and 39, respectively. The extension of the approach to meta-GGA functionals is described in this work.
The layout of this work is the following. Next, in section 2, we will outline the key pieces of FEM and present the optimal formalism for meta-GGA functionals with fractional occupations, reducing the problem into separate radial subproblems for each angular momentum l. Computational details including the studied selection of various LDA, GGA, and meta-GGA functionals are presented in section 3. Results of these functionals on the closed-shell and half-closed-shell atoms from H to Ar are presented in section 4. We will demonstrate that sub-μE h accurate total energies can be routinely determined with the new code for various well-behaved meta-GGA functionals, and that this allows the accurate determination of the truncation errors of various Gaussian basis sets. We will also show that taking full use of the symmetry inherent in the problem yields significant speedups in the calculations, enabling calculations to be performed at the CBS limit in a matter of seconds on commodity hardware. Finally, we will examine the density thresholds employed in the various density functionals considered in this work. The article concludes in a summary and brief discussion in section 5. Atomic units are used throughout, unless specified otherwise.

THEORY
In this section, we will give a brief overview of the theory necessary for implementing meta-GGAs with fractional occupations in nonrelativistic atomic calculations. We assume that the spin-σ orbitals are of the form where R σnl (r) are the spin-σ radial functions for primary quantum number n and angular quantum number l, and Y l m are complex-valued spherical harmonics. When each such spin− orbital is occupied by 0 ≤ f σnlm ≤ 1 electrons, the spin-σ electron density comes out as To achieve a spherically symmetric spin-σ electron density n σ (r) that only depends on the distance to the nucleus leads to reduction of eq 4 to the radial-only form of eq 5.
The radial orbitals are expanded in a numerical basis set χ μ (r) as and the energy is minimized in terms of the orbital coefficients C μn σ (l) . For brevity, we assume familiarity with the implementation of LDAs and GGAs in a finite basis set approach; 1 a detailed description of the procedure for atomic finite element calculations can be found in refs 34 and 37.
Substituting eq 7 into eq 5 leads to the compact expression n r P r r ( ) ( ) ( ) = (8) where the density matrix This mathematical structure was used in ref 34 to formulate an approach for LDA, GGA, and hybrid functionals (including Hartree−Fock theory) that reduces to solving a set of coupled radial eigenvalue equations, leading to significant savings in computational and storage requirements for the wave function.
In the following, we extend this approach to functionals that depend on the spin-σ local kinetic energy density τ σ and/or the density Laplacian ∇ 2 n σ as where ϵ xc is the energy density per particle that defines the used DFA and γ σσ′ = ∇n σ ·∇n σ′ is the reduced gradient. Note that DFAs (eq 11) are often also written in terms of the energy density f xc = nϵ xc as 2.1. Kinetic Energy Density Dependent Functionals. The positive-definite kinetic energy density τ σ , which is the most popular ingredient for meta-GGAs, is given by and Substituting eq 7 into eq 17 leads to our final expression Eq 18 is our first result: the local kinetic energy density can be rewritten as a sum of contributions from various angular momenta, which can be written solely in terms of radial density matrices. Furthermore, since we use the same radial basis set for all angular momenta, 37 eq 18 can be evaluated faster as where P σ was defined in eq 9 and we have introduced an angularweighted density matrix The final step needed for a SCF algorithm is the minimization of the total energy by variation of the orbital coefficients. 1 As all the DFA ingredients are now seen to be spherically symmetric, the density functional contribution to the energy from eq 12 simplifies to In the LCAO approach, varying the total energy E with respect to the orbital coefficients leads to the Roothaan equation F σ C σ = SC σ E σ , 1 where F μν σ = ∂E/∂P μν σ is the spin-σ Fock matrix, C σ and E σ are the orbital coefficients and the corresponding diagonal matrix of orbital energies, respectively, and S is the overlap matrix with elements S μν = ⟨χ μ |χ ν ⟩. 1 The use of spherical symmetry leads to the Roothaan equation splitting into radial subproblems F σ(l) C σ(l) = S (l) C σ(l) E σ(l) for every angular momentum l, where the radial Fock matrix is 22) and the radial overlap matrix is The LDA and GGA type contributions to the radial Fock matrix are independent of the angular momentum 34 which evaluates to which is the final piece of the implementation for τ-dependent functionals.

Density Laplacian Dependent Functionals.
The Laplacian of the density, ∇ 2 n σ , is straightforward to process, as the density is spherically symmetric by construction (eq 5). One must merely remember that the Laplacian has two terms in the spherical polar coordinate system: The Laplacian dependence leads to an additional contribution to the Fock matrix given by Note that the expression in eq 31, like the LDA and GGA contributions of eqs 24 and 25, applies to all angular momenta, at variance to the term arising from local kinetic energy density dependence in eq 27. Note also that the first two terms in the integral eq 31 can be combined with the GGA expression in eq 25, while the third term in eq 31 is of the same form as the first term in eq 27 and can be likewise evaluated together.

Numerical Stability
Close to the Origin. The local kinetic energy τ σ , eq 19, seems tricky to evaluate near the origin due to the second term of eq 15, that is, ∑ nl l(l + 1)n σnl (r)/r 2 . However, it is easy to see that the term is regular, as only s orbitals have electron density at the nucleus: the s orbital contribution is killed as l(l + 1) = 0, while n σnl (r) → 0 when r → 0 for l > 0.
A minor complication is that although n σnl (r) ≥ 0 by definition, in practice ∑ nl l(l + 1)n σnl can attain a small negative value due to finite numerical precision, which can be magnified by a large r −2 factor to generate a large negative contribution close to the nucleus. We have found that such cases only occur in the few quadrature points closest to the nucleus that carry small quadrature weights. For simplicity, we opted to stabilize this term simply by ensuring that it is non-negative by setting any negative contributions to zero.
The density Laplacian ∇ 2 n σ , eq 29, has a clear singularity at the origin in its first term, where the product of the numerically stable basis functions (see eq 35 below for discussion) and their derivatives is multiplied by a singular r −1 factor. The singularity is, however, integrable, as each term becomes regular at the origin when multiplied by r 2 of the quadrature weight. This also makes the Fock matrix expression of eq 31 regular, provided that ∂f xc /∂(∇ 2 n σ ) is finite when r → 0.

Finite Element Calculations.
All finite element calculations are performed with the HelFEM program, 34,37,44 which employs Libxc 35 to evaluate DFAs. The used HelFEM implementation is available in the public GitHub repository. 44 The numerical basis functions used in this work are defined in terms of piecewise polynomial finite element shape functions B μ (r) as The shape functions B μ in eq 32 are expressed within each element r ∈ [r i , r i+1 ] in terms of a primitive coordinate x ∈ [−1, 1] obtained with the transformation x r r r r r r r r ( ) 2 1, , As in our previous works on atomic calculations, 5,34,37,39 the shape functions are chosen to be Lagrange interpolating polynomials (LIPs) with the nonuniform nodes {x α } chosen from the Gauss− Lobatto quadrature rule, which avoids the Runge instability 45 and allows the use of very high-order numerical schemes. We have recently studied the use of Hermite interpolating polynomials (HIPs) instead of LIPs and found that HIPs and LIPs yield similar results with τ-dependent meta-GGA functionals. 38 Note that although eq 32 is numerically unstable for small r, we have recently shown that numerically stable basis functions are afforded by Taylor expansions of eq 32 for small r: 38 where R is a switching radius that can be chosen automatically by optimal matching of the left and right-hand sides of the piecewise definition in eq 35. We employ the numerically stable form given by eq 35 with high-order Taylor series matching the polynomial order of B μ (r) in all calculations of this work, as discussed in ref 38. Integrals are computed by Chebyshev quadrature with N quad points; a rule transformed to unit weight factor is employed for this purpose, as it provides nodes and weights in easily computable analytical form. 36,37,46 All calculations discussed in this work are converged with respect to the number of quadrature points.
The calculations employ 15-node LIPs (corresponding to 14th order polynomials), as this order was found to be sufficient for a rapid convergence of Hartree−Fock total energies in ref 37. The "exponential grid" of ref 37 is used in the present work with the default values r ∞ = 40a 0 and z = 2 for the practical infinity and the grid parameter, 37 as these values afford excellent accuracy in Hartree−Fock calculations. 37 We have recently shown that LDAs, GGAs, and meta-GGAs have similar grid requirements to those of Hartree−Fock. 38 The correctness of the present implementation has been verified by direct comparison with the general implementation described in refs 34 and 37. The energy and Fock matrix reproduced by the symmetry specialized version are in exact agreement with those from the general implementation.
The SCF calculations are started from orbitals obtained from a tabulated potential from a converged LDA exchange calculation, 2,34 and employ a combination of Pulay's direct inversion in the iterative subspace 47,48 (DIIS) and the augmented DIIS 49 (ADIIS) methods for reliable SCF convergence: we typically observe convergence within a dozen of SCF iterations in the fully numerical basis set.

Studied Atoms.
Although the presently used implementation can handle both heavy atoms and open shells (with the limitations of the presently examined nonrelativistic level of theory with a point nucleus and the use of fractional occupations), 34 for simplicity, we will examine the H, He, Li, N, Ne, Na, P, and Ar atoms at their ground-state configurations given in Table 1. With the exception of H and He, the same atoms were also considered in our recent study on the numerical  Table 1 that only feature fully occupied orbitals. This also enables the direct comparison of the present results to those obtained with molecular codes employing Gaussian basis sets, for example. The correctness of the present implementation is also obvious from the excellent level of agreement between the finite element calculations and ones performed with benchmarkquality Gaussian basis sets 39 that afford submicrohartree accuracy for light elements. The Gaussian-basis calculations in the AHGBS-9 39 and the aug-pc-4 basis 50−55 in fully uncontracted form (un-aug-pc-4) were performed with Erkale, 56,57 which likewise employs Libxc for DFA evaluation. All basis functions that are unnecessary to describe the ground states of Table 1 were removed from the Gaussian basis sets. In analogy to the finite element calculations, the Gaussian-basis calculations were started from error function fitted atomic LDA exchangeonly potentials. 5 3.4. Studied Density Functionals. Guided by the exploratory calculations and the work presented in ref 36, the density functionals considered in this work along with their literature references are shown in Table 2.
The meta-GGAs include the Tao−Perdew−Scuseria−Staroverov (TPSS) meta-GGA as well as its revision (revTPSS), both of whose correlation parts we have recently found to be numerically ill-behaved for alkali atoms at fixed electron density. 36 Next, the MS0, MVS, and SCAN functionals of Perdew and co-workers were included since they have been found to exhibit successively degenerating numerical behavior. 36 The SCAN functional has already been found to be illbehaved in fully numerical calculations by Bartoḱ and Yates, who suggested the rSCAN functional where the ill behavior is fixed by a well-behaved polynomial expansion. 72 rSCAN was then later used by Furness et al. 73 to build the r 2 SCAN functional. r 2 SCAN has showed extremely promising accuracy in applications; 90,91 however, it has been found to be ill-behaved in fully numerical calculations by Holzwarth et al., 25 who proposed another modification, the r 2 SCAN01 functional. (Note that we found the whole SCAN family to be ill-behaved in ref 36.) The TASKCC functional recommended by Lebeda et al. 92 is included as another recent meta-GGA, which appeared to be well-behaved in our recent studies. 36,38  All the meta-GGA functionals examined in this work depend only on τ. The reason for this is that Laplacian dependence is mainly included only in older meta-GGA functionals, which are well-known to be ill-behaved. Excluding such functionals leaves only functionals that depend on τ. Deorbitalized functionals, 93 which replace the τ dependence in recent DFAs with the density Laplacian through the use of a Laplacian-dependent kinetic energy functional are an exception to this rule; however, we have found many such functionals to be numerically ill-behaved already at fixed electron densities. 36 Moreover, it is well-known that kinetic energy functionals have singularities near the nucleus that pose issues for numerics and the stability of SCF calculations, 94 which are an issue even when pseudopotentials are used. 95 The basis function requirements for τ-dependent functionals were investigated in depth in ref 38, where it was found that the LIP basis reproduces the correct CBS limit even though it does not explicitly guarantee a continuous τ by construction unlike the Hermite interpolating polynomial (HIP) basis investigated in ref 38. We tentatively attributed this success to the action of the variational theorem: as discussed in ref 37, discontinuous derivatives would lead to higher kinetic energies, which are disincentivized by the variational minimization of the total energy.
Note that calculations on Laplacian dependent functionals� which are not considered in this work�should use at least a second-order HIP basis set in order to make ∇ 2 n well-defined everywhere. Such machinery was presented in ref 38, to which we refer for further details.

Convergence to the Basis Set and Density
Threshold Limit. We aim for total energies accurate to 0.1 μE h with respect to all parameteres in the present calculations. In practice, we determine that two calculations have converged to within 0.1 μE h precision if the energies E 1 and E 2 of the two calculations agree within 0.04 μE h : |E 1 − E 2 | < 4 × 10 −8 .
We therefore carry out a large number of calculations to determine reference total energies converged with respect to all the parameters controlling the calculation: in addition to the number of radial basis functions, we also study the convergence with respect to the number of quadrature points, as well as the effect of the density screening threshold employed in the DFA implementation.
When the values of the density functional and its derivatives are evaluated on the numerical quadrature grid by Libxc, 35 points with insignificant electron density n(r) ≤ ϵ as defined by a preset threshold ϵ > 0 are skipped, setting the energy density f xc and all its derivatives to zero. The rationale for such thresholding is grounded on the physical observation that both factors in the energy density f xc = nϵ xc (n, ...) decrease when n → 0 and therefore points with negligible density do not contribute meaningfully to the total energy; a practical issue is also that points with extremely small densities often yield divergent derivatives. Because of this, such thresholding is commonly used in density functional implementations.
This also extends further: typical density functional quadrature approaches discard basis functions with negligible values. For instance, the approach of Stratmann et al. 96 discards atomic basis functions χ α (r) in the sphere around the atom where r ( ) | | < with 10 10 = ; such an approach also leads to errors in the electron density, and the approach thereby relies on the small errors in the density not having a significant effect on the total energy.
A threshold ϵ = 10 −12 has been used in previous studies with HelFEM, which is also the default in the Psi4 97 and ORCA 98 programs, for example. In this work, we determine the density cutoff ϵ used in the evaluation of the density functional, specified with the Libxc function xc_func_set_dens_threshold, by considering a series of calculations performed with decreasing values in ϵ = 10 −n with n = 8, ..., 15. (The default value used in Libxc for most density functionals for three-dimensional systems is ϵ = 10 −15 , below which value numerical issues are often encountered due to the use of double precision arithmetic.) We deem convergence to have been reached at the largest value of ϵ for which the energy does not change any more.
We wish to underline here that the use of density thresholds in this work is not an approximation that is special to this work, as such finite thresholds are always used to ensure that the numerical implementation of any functional is stable. Instead, the density threshold should be considered part of the definition of the DFA, and a threshold that is small enough to reproduce the converged value should always be used. Unfortunately, the issue is that the values of the used thresholds are typically not discussed in articles suggesting novel DFAs. However, we do find that the studied functionals converge quickly in the density threshold.
We determine the convergence to the CBS limit by considering calculations with a sequence of increasing numerical basis sets composed of 5, 10, 15, 20, 25, 30, and 35 radial 15node LIP elements with with 69, 139, 209, 279, 349, 419, and 489 radial basis functions, respectively. Convergence to the CBS limit is established when the difference in total energy to the calculation with the largest number of radial basis functions is below the wanted precision. The convergence with respect to the quadrature is checked by doubling the used number of quadrature points, and checking whether the total energies of the SCF calculations performed with different quadratures agree.
All fully numerical calculations failed with the MVS and SCAN functionals, and these functionals were excluded from all analyses; the ill-behavedness of SCAN was already reported by Bartoḱ and Yates. 72 The obtained fully numerical reference energies for the remaining functionals are shown where available in Table 4 for H, He, Li, Be, and N and in Table 5 for Ne, Na, Mg, P, and Ar. Our analysis in Tables 4 and 5 also distinguishes cases where reliable reference energies were not achieved by the present computational scheme due to (i) failures with SCF convergence, (ii) failures with achieving the CBS limit, and (iii) floating point errors in the calculation.
As expected from our previous experience, Li, Be, Na, and Mg are a challenge for many functionals due to their loosely bound outer electrons, as can be seen from the large number of failed calculations for these atoms. Although the initial guess is good, as suggested by DIIS errors in the range of 10 −2 , we observe that many of the failed calculations go wrong already in the first iteration. For instance, the revM06, revM06-L, M06-SX, MN12-L, MN12-SX, MN15, B97M-V, and ωB97M-V calculations on Li start from a sensible total energy from the initial guess, but jump up thousands of Hartrees in energy upon the diagonalization of the first Fock matrix, which suggests that there are issues in these functionals' numerical behavior for some densities.
These observations motivate a systematical examination of the initial guess. The largest DIIS errors and largest energy changes in the first iteration that arise from the first full diagonalization of the Fock matrix are shown in Table 3. The functionals that stand out with a large DIIS error are MS0 and MVS, and all nine Minnesota functionals. The Ar atom has the largest initial DIIS error out of the studied atoms and functionals, with the exception of PW92 that encounters its largest initial DIIS error for the N atom.
Interestingly, the initial DIIS error does not appear to correlate strongly with the initial change in energy. Only B97, the Berkeley meta-GGA functionals B97M-noV and ωB97M-noV, and six out of nine Minnesota functionals show an increase in energy upon the first diagonalization; however, these increases are alarmingly large. We note again that the largest stability problems appear to be encountered with the Li and Na atoms, as can be seen from Table 3. Our baseline of HF, PW92, PBE, BLYP, B3LYP, and B97 converge without issues, with the exceptions of Li and N with B97. The B97 functional appears to be less smooth than the other baseline functionals, as evidenced by its need for more radial elements to reach the same convergence.
The TPSS, revTPSS and TASKCC functionals are wellbehaved, easily converging to the CBS limit for all atoms with the exception of Li, which required surprisingly many elements.
The MS0 functional is ill-behaved, failing to reach the CBS limit for Li, Be, Na, Mg, and P even with the extended numerical basis sets considered in this work.
The rSCAN functional is well-behaved, other than failing to reach the CBS limit for Li. In partial agreement with Holzwarth et al., 25 we find that the r 2 SCAN01 functional is better-behaved than r 2 SCAN, as calculations failed for Li and Na for the latter functional. However, we do not find evidence that r 2 SCAN01 is otherwise smoother, as the functional still takes more radial elements to converge than the well-behaved TPSS, revTPSS or TASKCC functionals, and as r 2 SCAN01 required more radial elements to converge P than r 2 SCAN did.
The ωB97X-noV range-separated hybrid GGA is wellbehaved and converges easily for all studied systems. In contrast, the B97M-noV and ωB97M-noV meta-GGAs fail to converge for Li, N, Na, and P. All of the failed calculations are characterized by large jumps in energy in the first iteration, as discussed above for several functionals with Li and Na.
The nine Minnesota functionals appear well-behaved for H, He, N, Ne, and Ar, as all nine functionals reach converged CBS limit energies. The M08-HX functional, however, requires many more radial elements than the other functionals, which can be explained by our recent observation in ref 38 of sharp nonphysical behavior in the functional. No Minnesota functional is successful for Li, while only revM11 is able to reach a converged CBS limit for Na. Only MN15, revM06, revM06-L, and M06-SX reach CBS limits for Be. For Mg, the CBS limit is reached by MN15, MN15-L, revM06-L, and revM11. All the studied Minnesota functionals except MN12-SX reach a CBS limit for P.
The converged density thresholds are shown in Table 6; the systems that failed to converge in Table 4 were excluded in this analysis. All functionals reach total energies converged to 0.1 μE h with a density threshold of ϵ = 10 −11 or larger, confirming the general reliability of the universally used screening approach, and confirming the reliability of the results previously obtained with the default threshold ϵ = 10 −12 .

Gaussian-Basis Truncation Errors.
Furnished with the converged fully numerical reference energies, we are able to determine truncation errors of Gaussian basis sets. This part of the study is partly motivated by our recent work in ref 99, where we observed unexpectedly large basis set truncation errors for hydrogen with the M06-L, 100 M11-L, 101 and B97M-noV functionals. We later found M06-L and M11-L to be illbehaved, 38 exhibiting large oscillations in the density Laplacian ∇ 2 n in the ground state of the hydrogen atom, which explains the large differences in energies in fully numerical and Gaussianbasis calculations.
The question of the accuracy of meta-GGA functional energies in Gaussian basis sets has not been addressed in the literature so far to the best of our knowledge. The recent study of Kraus 102 studied basis set extrapolations with modern density functionals using HelFEM, but does not appear to comment on the functional dependence of the accuracy in total energy.
We consider the aug-pc-4 basis set 50−55 in its fully uncontracted form (un-aug-pc-4), and our recent augmented hydrogenic Gaussian basis set (AHGBS-9). 39 The un-aug-pc-4 basis set has been optimized for the BLYP functional, 51,54,55 while the AHGBS-9 basis set and its polarized counterparts are constructed by considerations on one-electron ions, only. 39 The hydrogenic basis sets of ref 39 are large even-tempered basis sets aimed for benchmark accuracy calculations on atoms and molecules. Note that even-tempered basis sets are often used for studies on basis set completeness for their favorable properties. 103    The notation is analogous to that in Table 4.

Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article
The truncation errors for the un-aug-pc-4 basis set are shown in Table 7, and the errors for AHGBS-9 basis set are shown in Table 8. The Gaussian-basis energies were determined with the basis set truncation and density screening thresholds ϵ = 10 −12 and 10 12 = , respectively, and a (500, 974) quadrature grid. We observe that the truncation errors are strongly functional dependent, which is not surprising given the analogous FEM data in Table 4 and Table 5. Interestingly, even though un-augpc-4 has been optimized for the BLYP functional, it is often not the functional for which the smallest truncation error is observed: the lowest truncation error is often also observed for the PW92 or B3LYP functional. This is of course not surprising, as optimality of the exponents for a given functional does not prevent the truncation error for another functional with the same basis set being smaller.
In the case of the AHGBS-9 basis set, we observe that the smallest truncation errors are achieved for all atoms in case of HF calculations. This can likely be attributed to the quadratic character of the HF energy functional compared to the more complicated mathematical form of DFAs.

Examination into Differences
Between un-aug-pc-4 and AHGBS-9. AHGBS-9 generally yields smaller basis set truncation errors than un-aug-pc-4. This means AHGBS-9 is a better basis set, which is not surprising given its large size: the basis sets of ref 39 were designed for high-accuracy calculations on small systems.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article noV calculations on Li, and the TPSS, TASKCC, B97M-noV, ωB97M-noV, revM06-L, and M06-SX calculations on Be. The observed issues therefore only affect Li and Be, for which AHGBS-9 reproduces a lower total energy than un-aug-pc-4 in 16 calculations, and a higher one in 12.
In the former case, the largest decreases in total energy going from un-aug-pc-4 to AHGBS-9 are 15 μE h for Be with the MN15 functional, and 4.9 μE h for Li with the revTPSS functional. In the latter case, the largest decreases in total energy from AHGBS-9 to un-aug-pc-4 are 94 μE h for Li with the TASKCC functional, and 57 μE h for Be with the ωB97M-noV functional.
We shall investigate these discrepancies in the results further by examining the basis sets in detail. The examined configurations of Li and Be only have s electrons, and therefore examining the issue reduces to examining the s basis functions. Examining the s exponents, we see that they span the range 5.93 × 10 −3 a 0 −2 to 7.07 × 10 4 a 0 −2 in un-aug-pc-4 for Li, while the corresponding AHGBS-9 basis spans the range 7.05 × 10 −3 a 0 −2 to 6.85 × 10 5 a 0 −2 , suggesting that AHGBS-9 may not have sufficiently diffuse exponents in the case of Li. However, for Be we observe exponents in the range 1.11 × 10 −2 a 0 −2 to 1.39 × 10 5 a 0 −2 in un-aug-pc-4, while in AHGBS-9 the range is 6.51 × 10 −3 a 0 −2 to 1.22 × 10 6 a 0 −2 , which thus fully covers the range of exponents included in un-aug-pc-4 and suggests that the issue could be a lack of completeness within this range of exponents.
The completeness profiles 106 where α is a test function, μ and ν are atomic basis functions and ⟨μ|ν⟩ −1 denotes the μ, ν element of the inverse overlap matrix offer a visual tool to inspect the completeness of the studied basis sets. Gaussian functions with exponents α can be expanded exactly in the basis if Y(α) = 1, while functions that are orthogonal to the basis set have Y(α) = 0. The completeness profiles for the studied s functions in the Li and Be basis sets are shown in Figure 1. Figure 1 features oscillations in the profile of un-aug-pc-4 for large exponents α, while the profile for the large AHGBS-9 basis set is flat. These features are even clearer when examining the difference 1 − Y(α) in Figure 2: the AHGBS-9 basis appears considerably more flexible than un-aug-pc-4 in the case of Be in the same range of exponents.
We can thus summarize that a more complete Gaussian basis set may reproduce a higher energy than a smaller, less complete Gaussian basis set, depending on the functional. This indicates that the basis set truncation errors observed for various combinations of Gaussian basis sets and density functionals may have nontrivial dependence on the actual exponent values. The dependence of the optimal exponents on the DFA is caused by the differences in the optimal radial orbitals of various functionals. Our data in Tables 7 and 8 demonstrate that the use of a fixed Gaussian basis set can introduce functional dependent errors in total energies in the range of tens to hundreds of microhartree.
These findings again underline the importance of the present contribution in introducing fully numerical methods for the reliable determination of CBS limit energies. The issues discovered with the functional dependence of the exponents also underline another part of our discussion. For example, in the case of Be, using a converged PBE density to start a r 2 SCAN calculation shows that the PBE orbitals yield a total energy which is 0.957 mE h higher than that for the converged r 2 SCAN orbitals. As we have discussed in ref 12, NAOs are analogous to contracted basis sets, and this error is nothing but the contraction error made when using a minimal NAO basis for the PBE functional in a r 2 SCAN calculation. Although additional basis functions to allow breathing and polarization will allow for energy lowerings, the error of the minimal basis will reintroduce BSSE in calculations. We again underline that NAOs should be formed with the same DFA used in a polyatomic calculation to eliminate errors arising from differences in the optimal form of the radial orbitals. The resulting energies and timings are shown in Table 9. As HelFEM is a new project, the code has not been heavily optimized. Although many things could be done to optimize its performance, the code is usable in present form and fast enough to pursue investigations into the numerical stability of density functional approximations, for instance.
As the timings were obtained on the same machine with largely the same code, they give a good idea of the speedups achieved by the use of symmetry. The most time per iteration in the general program is spent on building the DFA components of the Fock matrix. This is also a part that experiences major speedups due to the use of symmetry to eliminate the angular degrees of freedom, as the quadrature over the solid angle is not needed, nor is the handling and pointwise evaluation of the spherical harmonics. The speedups are larger for DFT than for HF, and increase going from LDAs to GGAs to meta-GGAs.
Although obtaining results for heavy atoms with the general program may take up to tens of minutes, employing symmetry allows obtaining CBS limit results on commodity hardware in a matter of seconds.

SUMMARY AND DISCUSSION
We have presented the formalism necessary to implement meta-GGA functionals in atomic calculations within the finite element method, and implemented it in the free and open source HelFEM program. Furnished with the new implementation, we carried out a large number of calculations with 31 density functionals on the 10 closed-shell or half-closed-shell atoms from H to Ar to determine total energies converged to within 0.1 μE h with respect to all parameters controlling the calculation: the radial basis set, the quadrature scheme, as well as the density threshold in the density functionals' implementation in Libxc. 35 Excluding the nonconverging calculations, we found that a density screening threshold of 10 −11 a 0 −3 was able to reproduce total energies converged to 0.1 μE h for all studied functionals.
Ill behavior was observed in several density functionals. The Li and Na atoms proved to be the hardest systems in this study, which we attribute to their extended electronic structure. Pathological behavior was discussed for several functionals for the Li and Na atoms, where the diagonalization of a good initial guess results in thousand-hartree increases of the total energy. This points to issues with large derivatives, which were not examined in our recent study on the numerical behavior of density functionals, 36 and whose study was one of the central motivations of this work, as fully numerical calculations are stringent tests of density functionals' behavior.
Equipped with the fully numerical CBS limit energies, we proceeded to study basis set truncation errors in the AHGBS-9 39 and aug-pc-4 50−55 basis sets in fully uncontracted form (un-augpc-4). The truncation errors were found to be strongly dependent on the functional. Although AHGBS-9 is designed for benchmark studies and is thereby much larger than un-augpc-4, we found that un-aug-pc-4 afforded a lower total energy  than AHGBS-9 in 12 out of 28 calculations on Li and Be. (For all other systems, AHGBS-9 yielded systematically lower total energies.) Even though un-aug-pc-4 was found to have a more diffuse exponent than AHGBS-9 for Li, in the case of Be the unaug-pc-4 exponents were found to be included in the range of the exponents for AHGBS-9 and completeness profiles confirmed that AHGBS-9 is a more complete basis set. We therefore concluded that the use of fixed Gaussian exponents can introduce functional dependent errors in the range of tens to hundreds of microhartrees.
Our results underline the importance of fully numerical studies of novel density functionals. The timings presented in this work demonstrate that with the use of symmetry, functionals can be swiftly characterized by a fully numerical calculation. Furthermore, our implementation is open source and is freely available online for anyone for any purpose.
This study is the cornerstone on the road to employing modern finite element techniques for molecular calculations with NAOs. We hope to pursue along the path marked in ref 12 by introducing open source software for NAO calculations in upcoming work. However, as most algorithms required by a NAO program can be formulated independently of other technical choices made in the implementation, our plan is to pursue a modular approach. As