Robust Validation Of Approximate 1-Matrix Functionals With Few-Electron Harmonium Atoms

A simple comparison between the exact and approximate correlation components U of the electron-electron repulsion energy of several states of few-electron harmonium atoms with varying confinement strengths provides a superior validation tool for 1-matrix functionals. The robustness of this tool is clearly demonstrated in a survey of 14 known functionals, which reveals their substandard performance within different electron correlation regimes. Unlike spot-testing that employs dissociation curves of diatomic molecules or more extensive benchmarking against experimental atomization energies of molecules comprising one of standard sets, the present approach not only uncovers the flaws and patent failures of the functionals but, even more importantly, allows for pinpointing their root causes. Since the approximate values of U are computed at exact 1-densities, the testing requires minimal programming, and thus is particularly useful in quick screening of new functionals.


I. INTRODUCTION
The present implementations of the density matrix functional theory (DMFT), in which the one-electron reduced density matrix (a.k.a. 1-matrix) Γ(1 ′ , 1) plays the central role [1][2][3][4][5][6], are not on par in accuracy with the wavefunction-based methods of quantum chemistry and the parameterizations of the relevant functionals are far less sophisticated than those of their density functional theory (DFT) counterparts. Paradoxically, this state of affairs, coupled with the wealth of constraints that have to be satisfied by the 1-matrix functionals [7][8][9][10], bodes well for mathematical rigor and sound physical basis of approximate formulations of DMFT that are to emerge in the near future. Another favorable circumstance stems from the fact that, unlike in the case of DFT, the homogeneous electron gas does not provide a convenient starting point for construction and testing of DMFT-based approaches to the electron correlation problem. Consequently, more realistic model systems have to be employed in its place.
The two-electron harmonium atom, described by the nonrelativistic Hamiltonian [11,12] with N = 2, is an archetype of quasi-solvable systems of relevance to electronic structure theory. As such, it has been repeatedly used in calibration and benchmarking of approximate electron correlation methods, especially those involving DFT [13][14][15][16][17][18][19]. However, due to the simplicity of the expression for the ground-state energy of a two-electron system in terms of its 1matrix, the two-electron harmonium atom is of no interest to the developers of DMFT. In contrast, there are multiple reasons for which its congeners with N > 2 are ideal model systems in this context [20]. First of all, such atoms offer unlimited continuous tunability of extents and relative strengths of the dynamical and nondynamical electron correlation effects. Thus, for large values of ω, they are weakly correlated systems that, depending on N and the electronic state, are described by either one or a few Slater determinants [21]. Conversely, at the ω → 0 limit of a vanishing confinement strength, their electrons exhibit complete spatial localization and exclusively nondynamical correlation [11,12,[22][23][24][25]. The absence of a sudden Wigner crystallization allows for smooth interpolation between these weak-and strong-correlation regimes [26]. Second, lacking the electron-nucleus cusps in their wavefunctions, all harmonium atoms are amenable to calculations involving basis sets of explicitly correlated Gaussian functions that yield highly accurate electronic properties [27]. Indeed, energies of several electronic states of the three-and four-electron harmonium atoms are presently known within ca. 1 [µhartree] for arbitrary values of ω [28,29]. The respective 1-matrices and individual energy components are also available from such calculations. Third, exact asymptotics of these electronic properties are available at both the weak-and strong-correlation limits [21-25, 30, 31], facilitating imposition of well-defined constraints upon the DMFT expressions for the total energy. In summary, many-electron harmonium atoms are well suited for robust validation of approximate 1-matrix functionals. Although such a validation can be carried out in multiple ways, the most convenient (and computationally least expensive) approach relies on comparison between the exact and approximate correlation components of the electron-electron repulsion energy W . In Coulombic systems, partitioning of W parallels that of the 2-matrix [32,33]. Thus, W comprises the direct (Coulomb) and exchange components, given by the expressions and respectively, and the remainder U due to electron correlation. This reminder, which originates from the diagonal part of the 2-cumulant matrix, is the only contribution to the total energy not given by an explicit functional within DMFT. Consequently, evaluation of the difference between the (almost) exact U, obtained either from highly accurate calculations or asymptotic expressions, and its approximate counterpart computed with the exact Γ(1 ′ , 1) (originating from the same source as the exact U) offers a rigorous test of the performance of a given functional that involves minimal (if any) programming as no optimization of the 1-matrix is required.
In this paper, we employ the aforedescribed validation tool in a survey of the majority of the currently known 1-matrix functionals. The origins of their flaws and patent failures are elucidated.
The genuine "JKL-only" [34,52] functionals of the PNOF (Piris natural orbital functional) family obtain from model reconstructions of the 2-cumulant matrix and thus include the Coulomb integrals {J pq }, where J pq = ψ p (1)ψ q (2)|r −1 12 |ψ p (1)ψ q (2) , in the approximate formulae for U, alleviating the aforementioned deficiency. Unfortunately, the assumption of NOs possessing pairwise-identical spatial parts limits the applicability of the PNOF functionals to spin-unpolarized systems and species with all-parallel spins, each of those two cases requiring a distinct expression for U. For the sake of reader's convenience, these expressions are compiled below: 1. The PNOF1 functional [53,54]: for the spin-unpolarized systems and for the high-spin ones; 2. The PNOF2 functional [55]: (15) for the spin-unpolarized systems and for the high-spin ones; 3. The PNOF3 functional [56]: for the spin-unpolarized systems and for the high-spin ones; 4. The PNOF4 functional [57]: for the spin-unpolarized systems and for the high-spin ones; 5. The PNOF6 functional [58]: for the spin-unpolarized systems and for the high-spin ones, where In Eqs. (13)-(23), ν p equals one when the pth NO belongs to the set of the N natural spinorbitals with the greatest occupation numbers (where N is the number of electrons); otherwise it equals zero. The other quantities are defined as Two comments are in order here. First, because of its pairwise coupling of occupation numbers, the PNOF5 approach [59] is not included in the present study. Second, the presence of single-index terms in the expressions for U pertinent to spin-unpolarized systems implies the lack of invariance with respect to unitary transformations among NOs with degenerate occupation numbers -the same problem that afflicts some of the functionals of the first category.

III. DETAILS OF CALCULATIONS
The set of validation tools employed in the present work comprises the lowest-energy 2 P − and 4 P + states of the three-electron harmonium atom and the lowest-energy 1 D + , 3 P + , and 5 S − states of the four-electron species in combination with 19 magnitudes of the confinement strength ω that span the range from 10 −3 to 10 3 [28,29]. The three-electron wavefunctions of the 2 P − and 4 P + symmetries describe, respectively, the ground and the first excited states for all confinement strengths. On the other hand, the ground state of the four-electron harmonium atom has the 3 P + symmetry for large values of ω and 5 S − for small ones, the transition occurring at ω ≈ 0.0240919 [29]. The 1 D + state always lies above its 3 P + counterpart.

TABLES I-III HERE
In the case of the 2 P − and 4 P + states of the three-electron harmonium atom, the values of the correlation component U(ω) of the electron-electron repulsion energy accurate to within a few µhartree have been published elsewhere [31]. Application of the same numerical methods to the results of high-quality electronic structure calculations on the four-electron species [29] produces the data listed in Tables I-III. The computation of the NOs and their occupation numbers has been described previously [28,29].
At the ω → ∞ limit of vanishing electron correlation, the values of U(ω) tend to constants that equal [21,31] − 98 135 and for the 2 P − , 4 P + , 3 P + , and 5 S − states, respectively, reflecting the singlydeterminantal nature of the underlying wavefunctions. The fidelity with which these constants are reproduced reflects the performance of approximate 1-matrix functionals for systems of diverse spin multiplicities in which only the dynamical electron correlation is present.
In contrast, faithful reproduction of the analogous asymptotics for the . , indicates the suitability of a given functional for description of systems with the nondynamical electron correlation due to degeneracy of the zeroth-order wavefunction. Finally, the differences between the exact and computed values of U(ω) within the strong-correlation regime of small ω quantify the performance of approximate functionals for quasi-classical systems that lack strongly occupied NOs [22,25].
For the approximate expressions (4)- (12), the validation encompasses all the five states in question. On the other hand, for the functionals of the PNOF family, the test runs involve only the 4 P + , 1 D + , and 5 S − states due to the aforediscussed restrictions on the spin multiplicity.

IV. RESULTS AND DISCUSSION
Several dichotomies are discernible in the performances of the 1-matrix functionals for the 2 P − ground state of the three-electron harmonium atom (Table IV). Within the weak-correlation regime, the MBB, GU, BBC1, and BBC2 functionals reproduce the exact values of U with comparable accuracy. What sets them apart, however, is the ability to correctly describe the strongly correlated species. Here, the GU functional fares the best, whereas the BBC1 and BBC2 estimates for U exhibit a sudden drop in accuracy for ω < 10 −2 due to the "phase-switching" step functions present in Eqs. (8) and (9). The MBB, CA, and CGA expressions yield similarly poor approximations to U. Paradoxically, the performance of the ML and ML-SIC functionals, which is very unsatisfactory at the weak-correlation limit, improves significantly upon weakening of the confinement, making them the most accurate ones at the smallest values of ω.
Inspection of the data compiled in Table V reveals an entirely different situation for the 4 P + high-spin state where none of the functionals involving solely the exchange integrals is capable of producing reasonable approximations to U as ω → ∞. Their performance for the strongly correlated species is somewhat better with only the BBC1, BBC2, ML, and ML-SIC expressions being reasonably accurate.

TABLE VI HERE
Due to its multi-determinantal character, the 1 D + state of the fourelectron harmonium atom turns out to be particularly challenging for the "K-only" functionals. In fact, a rough agreement with the exact data is observed only for the GU, ML, and ML-SIC approximate values of U within the strong-correlation regime (Table VI). The positive-valuedness of the BBC1 and BBC2 estimates for small ω is worth noting in this context.

TABLES VII AND VIII HERE
According to the data displayed in Tables VII and VIII, the performance of the functionals under study for the 3 P + and 5 S − states of the four-electron species largely parallels that for the 2 P − and 4 P + states of its three-electron counterpart. Thus, within the weak-correlation regime of the 3 P + state, the MBB, GU, BBC1, and BBC2 functionals are again the most accurate (though somewhat less than in the analogous case of the 2 P − state). Similarly, the ML and ML-SIC functionals perform well only for very small values of ω. At that limit, they are more accurate then the GU approximation, whereas the MBB, CA, and CGA expressions again yield similarly poor estimates of U. The BBC1 and BBC2 functionals fare the worst, producing positive values of U.
Not surprisingly, the results of test calculations for the 5 S − high-spin state reveal extremely poor reproduction of the exact correlation components of the electron-electron repulsion energies by all the functionals, the only exception being the ML and ML-SIC ones at small confinement strengths. In contrast to the case of the 4 P + state, the BBC1 and BBC2 expressions yield positive values of U at the ω → 0 limit.

TABLES IX -XI HERE
The data compiled in Tables IX-XI uncover a substantial underestimation or even a complete neglect of the electron correlation effects in high-spin states by the functionals of the PNOF family. For the PNOF1, PNOF2, and PNOF3 ones, this flaw persists for the 1 D + state. On the other hand, among those included in the present survey, the PNOF4 and PNOF6 functionals are the only ones capable of describing the ω → ∞ limit of this state with decent accuracy. Although their performance quickly deteriorates upon weakening of the confinement, the latter functional yields estimates of U that are both negative and greater than the exact ones for all values of ω (Table X).
Overall, one finds the functionals under study disappointingly inaccurate. With the exception of the PNOF4 and PNOF6 ones, all of them fail for a system described by a wavefunction with two predominant Slater determinants. Moreover, none of the functionals performs satisfactorily for high-spin species within the weak-correlation regime, signaling an unbalanced description of systems with different spin multiplicities. At the strong-correlation limit, only the ML and ML-SIC expressions appear to perform reasonably well. However, it is not clear at this point whether their accuracy is not purely accidental in this instance.
Most of the flaws and failures uncovered by the present survey have readily traceable origins. The most obvious one is the presence of the "phaseswitching" step functions in the BBC1 and BBC2 expressions. Whereas seemingly improving description of dissociation limits of simple molecules with single bonds [41], it introduces discontinuities in the first-order derivatives of U(ω) with respect to ω and results in the quick deterioration of accuracy apparent upon further lowering of the confinement strength.
The poor performance of the ML and ML-SIC functionals for weakly correlated systems can be elucidated with equal ease. In such systems, U(ω) tends to a constant at the limit of ω → ∞ [see Eqs. (25)- (28)], whereas the electron-electron repulsion integrals grow like ω 1/2 . Consequently, the expressions premultiplying these integrals have to be asymptotically proportional to ω −1/2 . The deviations of the occupation numbers {n p } from their ω → ∞ limiting values (equal to either 0 or 1) scale asymptotically like ω −1 [12,60]. Hence, the leading asymptotic terms of the premultipliers have to be proportional to square roots of the products {n p n q } for the pairs comprising one weakly and one strongly occupied NO that exclusively contribute to U(ω) at this limit [61]. Obviously, the expressions (10) and (11) do not conform to this constraint.  ω) for U afforded by the MBB, GU, CA, GCA, BBC1, and BBC2 functionals, respectively. These asymptotics, which are compatible with the large-ω entries in Table VI, should be compared with the exact limit of 1 turns out to be given by the rather unwieldy expression ln(2 √ 454+10)−ln 33 ln 2 ≈ 0.672 996 that compares well with the actual data computed within the weakcorrelation regime (Table XII).
Finally, a comment on the uneven performance of the functionals for systems with diverse spin multiplicities is in order. This problem stems from the insufficiency of the occupation numbers alone for a proper reconstruction of the 2-cumulant matrix that is equally applicable to species with balanced and unbalanced spins [62].

V. CONCLUSIONS
A simple comparison between the exact and approximate correlation components U of the electron-electron repulsion energy of several states of fewelectron harmonium atoms with varying confinement strengths provides a stringent validation tool for 1-matrix functionals. The robustness of this tool is clearly demonstrated in a survey of 14 known functionals, which reveals their substandard performance within different electron correlation regimes. Unlike spot-testing that employs dissociation curves of diatomic molecules or more extensive benchmarking against experimental atomization energies of molecules comprising some standard set, the present approach not only uncovers the flaws and patent failures of the functionals but, even more importantly, also allows for pinpointing their root causes. Since the approximate values of U are computed at exact 1-densities, the testing requires minimal programming, and thus is particularly suitable for rapid screening of new functionals.
The conclusion that emerges from the survey of functionals reported in this paper is that the current approximate incarnations of DMFT are highly unlikely to compete with the more traditional approaches to the electron correlation problem. However, the present study clearly identifies the deficiencies that have to be rectified and provides obvious clues for the direction of the future work.