How good are recent density functionals for ground and excited states of one-electron systems?

Sun et al. [J. Chem. Phys. 144, 191101 (2016)] suggested that common density functional approximations (DFAs) should exhibit large energy errors for excited states as a necessary consequence of orbital nodality. Motivated by self-interaction corrected density functional calculations on many-electron systems, we continue their study with the exactly solvable $1s$, $2p$, and $3d$ states of 36 hydrogenic one-electron ions (H-Kr$^{35+}$) and demonstrate with self-consistent calculations that state-of-the-art DFAs indeed exhibit large errors for the $2p$ and $3d$ excited states. We consider 56 functionals at the local density approximation (LDA), generalized gradient approximation (GGA) as well as meta-GGA levels, also including several hybrid functionals like the recently proposed machine-learned DM21 local hybrid functional. The best non-hybrid functional for the $1s$ ground state is revTPSS. The $2p$ and $3d$ excited states are more difficult for DFAs as Sun et al. predicted, and LDA functionals turn out to yield the most systematic accuracy for these states amongst non-hybrid functionals. The best performance for the three states overall is observed with the BHandH global hybrid GGA functional, which contains 50% Hartree-Fock exchange and 50% LDA exchange. The performance of DM21 is found to be inconsistent, yielding good accuracy for some states and systems and poor accuracy for others. Based on these results, we recommend including a variety of one-electron cations in future training of machine-learned density functionals.


I. INTRODUCTION
Density-functional theory 1,2 (DFT) has become one of the workhorses of computational chemistry, material science, and related fields, as modern density functional approximations (DFAs) only require a reasonable amount of computational effort while providing a level of accuracy sufficient for semi-quantitative predictions on a broad range of systems. [3][4][5] Although hundreds of DFAs have been proposed, thus forming the infamous zoo of DFAs, 6 new DFAs continue to be developed in the aim to find more universally applicable DFAs that combine suitable levels of accuracy and numerical effort.
New DFAs can be constructed along various strategies. 5,[7][8][9] The traditional route to construct DFAs is to start from first principles and to impose known limits and constraints; this is the way along which many well-known functionals such as PBE 10 , TPSS 11,12 , and SCAN 13 have been constructed.
Semi-empirical fitting is another route for constructing DFAs. Here, the general idea is to introduce flexibility in the functional form by introducing several independent DFA components that are weighted by parameters which are optimized against some training dataset. Classical examples of semi-empirically fitted functionals include B3LYP, 14 Becke's 1997 functional 15 (B97) and several refinements thereof such as the HCTH [Hamprecht-Cohen-Tozer-Handy] functionals by Handy and coworkers, 16,17 a) Electronic mail: schwalbe@physik.tu-freiberg.de b) Electronic mail: kai.trepte1987@gmail.com c) Electronic mail: susi.lehtola@alumni.helsinki.fi as well as the Minnesota family of DFAs by Truhlar and coworkers 18,19 that has been reviewed by Mardirossian and Head-Gordon 20 .
In reality, the classification of functionals into ones built solely from first principles vs ones formed by semiempirical fitting is not always clear: for instance, the TPSS exchange functional is parametrized to yield the exact energy for the hydrogen atom's exact ground state density, 11,12 while the SCAN functional 13 includes parameters that are fit to data on noble gases. Modern semi-empirical functionals, [21][22][23][24][25][26] in turn, typically employ a combination of the two approaches by restricting the fits to known constraints.
DFAs from either route are widely used, given their suitable numerical accuracy and reasonable computational effort. However, the functionals obtained from the two routes tend to exhibit different behavior. For instance, while semi-empirical DFAs often deliver excellent descriptions of the total energy, they may fail to reproduce electronic densities of the same quality: a famous article of Medvedev et al. 27 initiated an intense debate about this in the literature, [27][28][29][30][31][32] ; it was even pointed out that any general mathematical measure of density error is too arbitrary to be universally useful. 33 DFAs built on physical first principles, in contrast, often yield steady performance in a variety of applications, but may not achieve the same level of accuracy as tailored functionals for specific types of systems.
One of the most important limitations of present-day DFAs, regardless of their design, is the self-interaction error (SIE): an artificial interaction of the electrons with themselves. This error is related to density delocalization error and the fractional electron problem, 34,35 and leads to incorrect dissociation limits 36 and barrier heights, 37 for instance. Recent avenues for circumventing SIE in DFAs involve determining the electron density with another method, such as Hartree-Fock 38,39 or multiconfigurational wave function theory 40,41 . Other types of approaches have also been proposed in the literature. To solve the self-interaction problem, Perdew and Zunger 42 (PZ) proposed an orbital-by-orbital self-interaction correction (SIC) where E KS is the Kohn-Sham (KS) energy functional 2 and the self-interaction error (SIE) is defined by Here, n iσ is the electron density of the i-th occupied orbital with spin σ and E J and E xc denote the Coulomb and exchange-correlation energy functionals, respectively. The idea behind PZ-SIC is that the selfinteraction error defined by equation (2) vanishes for the exact functional, 42 and thereby the Perdew-Zunger functional of equation (1) is a better estimate for the total energy than the uncorrected Kohn-Sham DFA E KS ; indeed, the PZ functional is exact for one-electron systems such as the H + 2 molecule with approximate DFAs. Despite the simple logic used to construct the PZ-SIC functional, the PZ-SIC method turns out to be quite complicated. The introduction of the explicit orbital dependence in equations (1) and (2) breaks the unitary invariance of the energy functional, 43 requiring costly unitary optimization of the orbitals (see Ref. 44 for discussion). However, even though the resulting method is known to correct charge transfer errors and barrier heights, it does not lead to improved atomization energies with GGA and meta-GGA functionals in general. 45 Continued research has illuminated other important theoretical aspects of PZ-SIC. First, the orbital dependence in equations (1) and (2) has been recently shown to require the use of complex-valued orbitals for proper minimization, as real-valued orbitals can be shown to correspond to high-order saddle points. 46 When complexvalued orbitals are employed, the total energy is lowered, and PZ-SIC does lead to improved atomization energies for some GGA functionals; however, more accurate atomization energies can be obtained at significantly smaller cost with several standard DFAs. 47 Second, the orbital dependence in equations (1) and (2) has also been shown to lead to the existence of several local minima in the orbital space. 46 This problem has been recently shown to persist also in a related SIC method 48 based on the use of Fermi-Löwdin orbitals (PZFLO-SIC), where various choices for the orbital descriptors lead to distinct local electronic minima. 49 The existence of such local minima is a significant and underappreciated aspect of PZ-SIC and PZFLO-SIC calculations, as finding the true ground state may require extensive sampling of the space of the various possible localized electronic configurations or bonding situations.
Despite their theoretical shortcomings, PZ-SIC and PZFLO-SIC have been found useful in many applications [50][51][52] and we are positive that several of the aforementioned issues in PZ-SIC and PZFLO-SIC can be addressed by developments in the related theories by changing the way the self-interaction correction is applied. One possible way to achieve improved results would be to revisit DFAs based on the requirements of SIC calculations. 53 It is known that present-day DFAs yield poor estimates for the noded electron densities that are involved in SIC calculations. 54,55 Sun et al. 54 demonstrated that the ground and excited state densities of the hydrogen atom (as well as of H + 2 , see below) lead to large relative errors in the exchange-correlation energy compared to the exact values, but we are not aware of any self-consistent calculations on this issue.
Following recent discussion in the literature on the accuracy of DFAs on the electron densities of small atoms and ions [27][28][29][30][31][32][33] and motivated by the obvious connection of one-electron errors (OEEs) to the PZ-SIC and PZFLO-SIC methods, in this work we will analyze the OEE of various functionals for the 1s ground state as well as the 2p and 3d excited states of hydrogenic ions Z (Z−1)+ , whose exact energies are well-known to be given in atomic units by where Z is the atomic number, n ≥ l + 1 is the principal quantum number, and l is the angular momentum. As was mentioned above, calculations of ground and excited states of the hydrogen atom and of the 1σ g ground state and 1σ u excited state of H + 2 have been discussed by Sun et al. 54 with non-self-consistent electron densities, while the 1s ground states of hydrogenic mononuclear cations as well as the 1σ ground states of hydrogenic diatomic cations have been discussed recently by Lonsdale and Goerigk 56 using self-consistent calculations. The novel contribution of this work is to address (highly) excited states with noded electron densities of hydrogenic cations self-consistently. Importantly, like the 1s ground state, the 2p and 3d excited states (as well as the analogous higher excited states like 4f ) are the lowest states of the corresponding symmetry, and the groundstate Kohn-Sham scheme is applicable to such excited states as well as shown by Gunnarsson and Lundqvist 57 .
We pursue thorough density functional investigations of the 1s, 2p, and 3d states of hydrogenic ions in benchmark-quality Gaussian basis sets specially suited for this purpose with a selection of 56 popular DFAs, including the recently developed, highly sophisticated machine-learned DeepMind 21 (DM21) local hybrid functional. 58 The layout of this work is as follows. The computational details are presented in section §II, and the results are given in section §III. A summary of our findings and an outlook for further investigations is given in section §IV. Atomic units are used throughout, unless specified otherwise.

II. COMPUTATIONAL DETAILS
We only use free and open-source software (FOSS) in this work, following the philosophy discussed in Ref. 59. PySCF 60 is an electronic structure code for all-electron calculations using Gaussian-type orbitals (GTOs). As we are targeting one-electron states of specific symmetry (s, p, or d states), following Gunnarsson and Lundqvist 57 we truncate the basis set in all calculations to contain functions only of the pursued symmetry: calculations on the 1s/2p/3d state only include the basis functions of the corresponding symmetry (s, p, or d functions, respectively) from the chosen parent basis set. This procedure has two important features: the 2p and 3d excited states become the ground state in the reduced-basis calculation, and the computational requirements are smaller since fewer integrals need to be calculated in the reduced basis than in the original basis set.
The one-electron guess-which is exact for oneelectron systems and thereby is also expected to be accurate for calculations employing DFAs as well-is used in all calculations. 61 To ensure that the SCF procedure converges to the global minimum instead of a saddle point, the following procedure was used. First, a regular SCF calculation was performed with PySCF with default settings; direct inversion in the iterative subspace (DIIS) is used to accelerate these calculations. 62,63 Next, convergence to saddle point solutions was checked: cases where the SCF converged to a final energy higher than that of the initial guess were restarted, with new calculations employing iterative diagonalization with level shifting 64 instead of DIIS to converge to the ground state. All calculations reported in this work are fully converged to a threshold of 1 × 10 −7 E h .
For the GTO basis sets, we use the family of hydrogenic Gaussian basis sets 65 (HGBS-n) that have been designed for high-accuracy calculations on atoms and small molecules. A special feature of the HGBS basis sets is that the basis for atomic number Z is determined by a universal even-tempered basis set for the ions Y (Y −1)+ for Y ∈ [1, Z], whereas augmented hydrogenic Gaussian basis sets (AHGBS-n) add further functions for describing the Z = 0.5 one-electron ion. 65 The parameter n controls the relative precision of the hydrogenic Gaussian basis; (A)HGBS-n reproduces the exact total energies of the one-electron ions to an approximate relative accuracy of 10 −n . 65 The motivation of this approach in Ref. 65 was that a many-electron atom experiences a screened nuclear charge that can be rewritten in terms of a radially dependent effective charge Z eff = Z eff (r) with the asymptotic limits Z eff (0) = Z and either Z eff (∞) = Z ∞ with the asymptotic limit Z ∞ = 0 for Hartree-Fock and DFT or Z ∞ = 1 for the exact effective potential. 61 Another feature of the HGBS basis sets is that the functions of various angular symmetries are determined independently of each other, which facilitates the formation of polarized counterparts of the basis sets that are essential for studying molecules and excited states, as additional shells are added to the basis like lego blocks. 65 Following Ref. 65, the basis set with p ≥ 1 polarization shells and accuracy n is denoted (A)HGBSPp-n. The definition of polarization shells varies by atom (see Ref. 65 for discussion); however, as we only include the functions of the pursued symmetry in each calculation, the choice of the polarization level of the (A)HGBSPp-n basis set does not matter as long as the original basis contains functions of the highest targeted angular momentum for the targeted atom, that is, d functions in this work. For the reasons listed above, the hydrogenic Gaussian basis sets of Ref. 65 are ideally suited for the present study-as will be demonstrated in section §III by benchmarks with functions from the polarization consistent (pc-n) basis sets 66 and their augmented versions 67 (augpc-n)-and, as will be discussed in section §III A, we will take the exponents from the AHGBSP3-n basis sets in this work. All basis sets were taken from the Basis Set Exchange. 68 The Libxc library 9 -which implements over 600 DFAs-is used in PySCF to evaluate the DFAs. The library provides access to a vast variety of DFAs, of which 49 were chosen for this work; see table 1 for the complete list of investigated functionals. Our selection includes functionals of the first to the fourth rung of Jacob's ladder, 123 that is, local density approximations (LDAs), generalized gradient approximations (GGAs), meta-GGAs, as well as global and range-separated hybrid functionals. In addition, we consider six hybrids of rSCAN and r 2 SCAN with varying fractions of Hartree-Fock exchange discussed in Ref. 120; these functionals were defined in the PySCF input files as weighted combinations of r( 2 )SCAN exchange and Hartree-Fock exchange + 100% r( 2 )SCAN correlation. The DM21 functional was also chosen for this study; we use the original implementation in PySCF of Kirkpatrick et al. 58 . This brings up the total to 56 functionals for this study. An unpruned (300,590) quadrature grid is used in all calculations, including the non-local correlation component in B97M-V, ωB97M-V and LC-VV10.
As the total energies scale as E n ∝ Z 2 according to equation (3), the results will be analyzed in terms of absolute relative errors (AREs). Hydrogenic estimates show that the approximated exchange-correlation energy scales like Z in the large Z limit, 124 meaning that the relative errors should tend to zero like 1/Z. The ARE for a given state of a given ion is given by The information in the AREs is analyzed with two further error metrics. The mean state error (MSE) measures the overall functional error over all ions by averaging the   ARE over all ions

c -2 p c -3 p c -4 a u g -p c -2 a u g -p c -3 a u g -p c -4 u n c -a u g -p c -2 u n c -a u g -p c -3 u n c -a u g -p c -4 A H G B S P 3 -5 A H G B S P 3 -7 A H G B S
The overall error (OE) for a functional is obtained by further averaging the MSE over all considered states (1s, 2p, and 3d) III. RESULTS

A. Basis set convergence
Before pursuing density functional calculations, we analyze the basis set truncation errors (BSTEs) for the oneelectron cations in the polarization consistent and hydrogenic Gaussian basis sets. We aim for a mean BSTE smaller than 5 × 10 −5 E h for the whole benchmark set ranging from H 0 to Kr 35+ to ensure that our results are converged close to the complete basis set limit.
Unrestricted Hartree-Fock (UHF) is exact for oneelectron systems and thereby gives the exact energy E UHF n in the studied basis; the difference of E UHF/basis n and the exact analytical energy (equation (3)) is therefore a variational measure of the BSTE for the state with given n of the studied hydrogenic ions. The calculated mean BSTEs for a variety of polarization consistent and hydrogenic Gaussian basis sets are shown in figure 1; additional results can be found in the supplementary material. Unsurprisingly, uncontracting the (aug-)pc-n basis sets-yielding the unc-(aug-)pc-n basis sets-results in a noticeable decrease of the BSTE, because the contractions were determined in Ref. 66 with the BLYP functional 70-72 that suffers from SIE for the 1s state, while the p and d functions in the basis set describe either polarization effects or the occupied p or d orbitals in the screened neutral atom. Although the large uncontracted polarization consistent basis sets exhibit satisfactory performance for the 1s state, they result in much larger errors for the 2p and 3d states; this error is again caused by the p and d orbitals in the neutral atom being screened by the core electrons, which results in the lack of tight p and d basis functions that are necessary for the 2p and 3d states of the one-electron ions.
In contrast, the primitive (not contracted) hydrogenic Gaussian basis sets of Ref. 65 show uniform accuracy for the 1s, 2p, and 3d states, and as can be observed in figure 1, the targeted mean BSTE threshold is roughly achieved already with the AHGBSP3-7 basis set. The AHGBSP3-9 basis sets yield errors below the desired threshold for all states, and is therefore chosen for all the remaining calculations of this study.
Although this analysis was limited to Hartree-Fock calculations, we note that the basis set requirements of Hartree-Fock and DFT are known to be similar. 125 Furthermore, reliable reference energies for DFAs can be obtained with fully numerical methods, [126][127][128] and exploratory calculations presented in the supplementary material confirm that the BSTEs in the AHGBSP3-9 basis are small also for DFAs.

Exploratory analysis
We begin the analysis by a graphical study of the results of the SPW92, PBEsol, revTPSS, MN15-L, BHandH and DM21 functionals in figure 2. This limited set of functionals contains LDA, GGA, and meta-GGA functionals from first principles (SPW92, PBEsol, and revTPSS, respectively), semiempirical functionals (MN15-L and DM21) as well as hybrid functionals (BHandH and DM21).
As will be discussed in section §III B 2, revTPSS is the most accurate meta-GGA functional for the 1s state. In figure 2, revTPSS is outperformed by DM21 only for He + , and otherwise revTPSS affords much lower errors than the five other functionals in the figure. In contrast, the performance of DM21 is inconsistent. DM21 has lower errors for light ions than for heavy ions, but the curve is kinked for the light ions. DM21's curve becomes smooth for heavy ions, but DM21 is also less accurate for heavy ions. MN15-L also shows a kinky behavior with lower errors for light ions; these non-systematic features of DM21 and MN15-L can be tentatively explained by their semiempirical character; the curves for the first principles functionals are smoother.
The functional errors for the 2p state are shown in figure 2(b). The performance for the 2p state is strikingly different compared to the 1s state shown in figure 2(a). The plots for the 2p state in figure 2(b) show more structure and curve crossings. The behavior of  DM21 is qualitatively different from that of the other functionals: DM21 shows large relative errors for light atoms and lower relative errors for heavy atoms, while most of the other functionals shown behave similarly to each other. The only other exceptions to this are the SPW92 and MN15-L functionals that show dips at Z 3 and Z 4, respectively; the two functionals are thus oddly more accurate for some values of Z than others.
The errors for the 3d state are shown in figure 2(c). BHandH has small errors for all ions for the 3d state. The behavior of DM21 and MN15-L again differs qualitatively from the other functionals. While DM21 shows less variation for the 3d state than for the 2p state, MN15-L does the opposite: MN15-L has large errors for light ions, becomes nearly as accurate as BHandH for Z 22, while the relative error increases again for heavier ions.

Full analysis
The MSEs and OE for all studied functionals are shown in  the 1s, 2p and 3d states as well as in terms of the overall error can be found in the supplementary material. Clearly, the performance of all LDAs is practically identical. This suggests that the functional error for LDAs is limited by the simple functional form. While LDAs show larger errors than GGAs and meta-GGAs for the 1s state, they perform better than GGAs and meta-GGAs for the 2p and 3d states.
PBEsol 104 is the GGA that yields the smallest errors for the 2p and 3d states, as well as the smallest overall error. The XLYP GGA has a lower error than PBEsol for the 1s state. Although XLYP and even PBEsol are better for 1s states than any LDA, they have higher OEs than any LDA because of their considerably poorer performance for the 2p and 3d states. Analogous findings apply also to all other studied GGAs.
The best meta-GGA for the 1s state is revTPSS, 99,100 closely followed by rSCAN, 121 and r 2 SCAN 118 (see table 1 or the supplementary material). The best meta-GGA in terms of overall error is MN15-L. 97 Hybrid functionals have better accuracy, as they contain some Hartree-Fock exchange which is free of selfinteraction. The best hybrid GGA functionals in terms of overall error are BHandH 73 and QTP17. 109 BHandH has low MSEs for all states and has the best overall performance, which can be understood by its composition of 50% of Hartree-Fock exchange and 50% LDA exchange + 100% Lee-Yang-Parr correlation. QTP17 has the second best performance for all states; it, too, contains a mixture of Hartree-Fock (62%) and LDA exchange (38%).
The best functionals of each rung in terms of overall error are SPW92, PBEsol, MN15-L, and BHandH, respectively. The corresponding error distributions are summarized in comparison to DM21 in figure 3. Interestingly, the performance of the DM21 functional appears similar to that of MN15-L.
Following Medvedev et al. 27 , the calculated OE for all functionals and ions is plotted against the publication year in figure 4. As is clear from this plot, the improvement in one-electron error is not fully systematic and features a significant amount of spread and some no table outliers   Our trends and absolute values for the MSE for the 1s state are in satisfactory agreement for the subset of functionals studied in both works, although we did identify basis set incompleteness issues in some results of Ref. 56 as discussed in the supplementary material. The largest basis set incompleteness effects are observed for the M11-L and M06-L Minnesota functionals, which are known to converge remarkably slowly to the basis set limit. 131 Lonsdale and Goerigk 56 only studied one LDA functional (SVWN); we considered more LDAs and found them to have similar performance. Lonsdale and Goerigk 56 found OLYP to be the best GGA functional for the 1s state; we also considered XLYP and found it to yield a considerably lower MSE for the 1s state than OLYP. Lonsdale and Goerigk 56 included a broader set of hybrid functionals separating global, range-separated hybrids and double hybrids; however, our main motivation is the connection to self-interaction corrected methods where hybrid functionals are typically not used. We found rSCAN50/r 2 SCAN50 to be the best hybrid functional for the 1s state, while Lonsdale and Go-erigk 56 determined TPSSh and SCAN0 to be the best hybrids. All rSCAN and r 2 SCAN based hybrid functionals, i.e., rSCANh/r 2 SCANh, rSCAN0/r 2 SCAN0, and rSCAN50/r 2 SCAN50 as well as TPSSh have a good performance for the 1s state. The revTPSS functional is found in our work as well by Lonsdale and Goerigk 56 to be the best non-hybrid meta-GGA functional for the 1s state.

IV. SUMMARY AND DISCUSSION
We used exactly solveable hydrogenic cations in their 1s ground state and 2p and 3d excited states to determine the self-consistent one-electron error for 56 density functionals including the novel DM21 of Kirkpatrick et al. 58 , employing the methodology of Gunnarsson and Lundqvist 57 for the excited state calculations. In accordance with an earlier finding by Sun et al. 54 apparently based on non-self-consistent calculations for the hydrogen atom and molecule and one LDA functional, we find for 36 hydrogenic cations that all LDAs perform better for the excited 2p and 3d states than any of the tested GGAs and meta-GGAs. The performance of various LDAs appears to be almost identical, as the calculated errors are nearly indistinguishable, suggesting that the errors are limited by the simple functional form used in LDAs. Sun et al. 54 pointed out that larger errors for excited states are a necessary consequence of orbital nodality.
The revTPSS functional is the best performing meta-GGA for the 1s state, tightly followed by the rSCAN and r 2 SCAN functionals. MN15-L shows a better overall performance than LDAs for all states; however, the performance of MN15-L is non-systematic like that of DM21.
Hybrid functionals like BHandH and QTP17 have the best overall performance as they explictly include some fraction of Hartree-Fock exchange. Moreover, both BHandH and QTP17 are mixtures of Hartree-Fock and LDA exchange, leading to the good observed accuracy.
DM21 turns out to be only close to exact for the 1s state OEE from H 0 to B 4+ (see figure 2). For He + to B 4+ DM21 shows also good performance for 2p and 3d states. However, over the whole range of investigated species H 0 to Kr 35+ , DM21 exhibits various trend changes and an overall inconsistent performance. Thus, one might improve the next generation of the DM21 functional by including more one-electron cations in the training sets for various elements in the periodic table. This might increase the consistency of promising properties of such kind of machine-learned functionals.
We found PBEsol to be the most accurate GGA functional for the 2p and 3d states. PBEsol is also the most accurate GGA functional overall. These findings are interesting to contrast with that of Lehtola, Jónsson, and Jónsson 47 , who showed that PBEsol is one of the few functionals whose accuracy improves when PZ-SIC is applied with complex orbitals. The development of novel DFAs with reduced one-electron error could therefore be useful for PZ-SIC calculations, as the reduced oneelectron errors (equation (2)) would affect the numerics of the PZ correction (equation (1)) and might alleviate wellknown issues with PZ-SIC and PZFLO-SIC discussed in section §I.
Note added in proof After the acceptance of this paper, we became aware of a preprint by Lonsdale and Goerigk 132 that includes discussion on excited states of hydrogenic cations.

SUPPLEMENTARY MATERIAL
Exploratory finite element studies of basis set truncation errors in density functional calculations. Sorted rankings of the functionals by errors for the 1s, 2p and 3d states as well as the overall error. Bar plots of the errors for all studied functionals. Comparison of the 1s data to the study of Lonsdale and Goerigk 56 with additional basis set incompleteness studies. Tables of the calculated total energies for the 1s, 2p, and 3d states for all studied functionals.