Non-empirical, low-cost recovery of exact conditions with model-Hamiltonian inspired expressions in jmDFT Free

Approximate density functional theory (DFT) is one of the most widely used electronic structure methods throughout chemistry and materials science.^1–3 Presently available pure exchange-correlation (xc) approximations in Kohn-Sham DFT suffer from both one- and many-electron self-interaction errors (SIEs),^4–8 more commonly referred to as delocalization error (DE).^9–11 The manifestations of the delocalization error^12–14 include eroded accuracy for key properties for which DFT is often used, such as dissociation energies,^4–8 barrier heights,^5,15–18 band gaps,¹⁹ electron affinities,^20,21 and spin-state ordering.^22–28 Despite these shortcomings, the favorable balance of computational cost and accuracy that makes DFT tractable for the optimization and dynamics of larger (i.e., up to thousands of atoms) molecules and materials motivates its continued improvement.

A large number of strategies have been pursued^29–35 to improve DFT energetics broadly. Another class of approaches has explicitly focused on DE elimination in especially SIE-plagued systems (e.g., d- and f-valence materials) through explicit removal of the Kohn-Sham orbital self-interaction error^36–41 or other schemes to localize orbitals.^9,42,43 Generalizations⁴⁴ to Kohn-Sham DFT, e.g., global or range-separated hybrids^45–55 and DFT+U,^56–60 fall between these two limits. The efficacy of these approaches in improving materials and molecular properties has been explained^{4,7,11,20,59–70} by noting that they generally reduce DE. A number of variants^71–73 of locally focused DFT+U-like methods have been developed in recent years. A general shortcoming of most of these approaches is that although DE may be eliminated, the simultaneous mitigation of other errors, such as density^74,75 or static correlation error (SCE),⁷⁶ has eluded all but a few methods.^35,71,72,77 In the present work, we focus on the relationship between these errors and extend upon the DFT+U-inspired judiciously modified DFT (jmDFT) approach⁷¹ to simultaneously correct static correlation and delocalization errors. A limitation of the initial jmDFT effort⁷¹ was that coefficients were obtained by fitting to the errors of the functional with respect to an exact flat plane. In this work, we eliminate that remaining empiricism by introducing fully nonempirical coefficients. We also demonstrate the extension of the nonempirical approach to higher angular momentum.

We recently introduced⁷¹ the jmDFT approach in which (i) a specific condition of an exact electronic structure method is defined, (ii) deviations of an exchange-correlation functional (e.g., semilocal) in DFT from this exact behavior are obtained, (iii) corrections are fit to low order functions of the error in (ii), and (iv) those corrections are incorporated into the self-consistent DFT energy evaluation. In both this and prior work,⁷¹ we employ the flat plane^78–80 condition as the feature of an exact electronic structure method we aim to recover from semilocal DFT.

The flat plane is the union of two constraints: piecewise linearity with addition of fractional charge^56–59 and unchanged energy^81,82 when an electron is shifted from a spin up (α) to an isoenergetic spin down (β) orbital with constant charge (Fig. 1). Behavior along the flat plane has also been used to assess wavefunction theory methods^83,84 and in functional development.^77,85 As in prior work, we define the cationic fractional charge line (FCL+) as the line that connects the N-1- and N-electron configurations, e.g., by populating an n_α orbital, and the neutral FCL0 as the line that connects the N- and N+1-electron configurations, e.g., by populating an n_β orbital (Fig. 1). The FCL+ and FCL0 should be linear in fractional electron count, and their slopes correspond to the ionization potential (IP) and electron affinity (EA), respectively, of the N-electron system. This difference in IP and EA gives rise to a derivative discontinuity^79,80,86,87 at the point where the two FCLs meet. Over all pairs of (n_α, n_β), the flat plane connects numerous FCLs, and the isoenergetic line of the N-electron system is referred to as the fractional spin line (FSL, see Fig. 1).

Analysis of experimental ionization potentials of gas phase ions⁸⁸ revealed⁸⁹ that nearly all electron configurations give rise to this flat plane shape, but special combinations of excited states occasionally result in different shapes that lack the seam at the FSL, i.e., by having comparable IP and EA. As in prior work,⁷¹ we focus on electron configurations in which an initially empty valence shell in the N-1-electron ion (e.g., singlet s⁰) is populated first by an n_α electron to form the N-electron ion (e.g., doublet s¹) followed by an n_β electron to form the N+1-electron ion (e.g., singlet s²). For some of the cases studied here, this transition corresponds to a neutral N-electron ion (i.e., H, Li, Na) and N+1-electron anion, and for others (i.e., He⁺, Be⁺, Mg⁺), it corresponds to a neutral N+1-electron configuration.

In the present work, we extend beyond the s shell valence systems that have been studied in detail^71,82,85 to p and d electron configurations. To define a fractional spin error in these systems, half-filled p or d shells have been studied previously by simultaneously shifting all orbital occupations from spin up to spin down and computing changes in energy with fractional spin.^81,82 Due to our interest in developing nonempirical coefficients for jmDFT, we instead study analogues to the s valence shell ions, i.e., N-1-electron singlet p⁰ or d⁰, N-electron doublet p¹ or d¹, and N+1-electron singlet p² or d². The singlet N+1-electron configurations are no longer ground states as they were in the s valence cases but are well defined⁹⁰ as the lowest singlet state for each ion. The N-electron references here are selected to ensure that only p and d orbitals are emptied and filled in the flat plane endpoint configurations, i.e., N-electron states of B, C⁺, and N²⁺ for p¹ ions and Cr⁵⁺ for d¹ ions.

Evaluation of the flat plane error in an exchange-correlation (xc) functional (e.g., semilocal, generalized-gradient approximation or GGA) motivates possible functional forms for model-Hamiltonian-inspired corrections (Fig. 2). We previously showed that over a range of 23 ions and numerous functionals (e.g., local density approximations, GGAs, and hybrids), the shape of the flat plane error from DFT is generally similar (Fig. 2). Regions of concavity near the FCLs are joined by a region of convexity at the FSL (Fig. 2). In comparison with GGAs, there is more FSL convexity for hybrids and DFT+U. The nature of the error at the N-electron border suggested⁷¹ that low order corrections with few coefficients could be used to correct the majority of the flat plane error. These observations led us to compare⁷¹ five functional forms that could eliminate the flat plane error in jmDFT, three of which were sufficiently flexible to eliminate the flat plane error across a range of s valence ions and molecules.

A shortcoming of the initial jmDFT implementation was the need to carry out calculations at fractional numbers of electrons, which is not possible in many electronic structure codes. Here, we develop nonempirical coefficients that can be obtained from integer electron properties. To obtain coefficients from first-principles, we focus on the subset of jmDFT functionals in which semilocal DFT errors at integer electron count are not corrected.⁷¹ These forms have close connection to the widely applied, rotationally invariant and simplified form of DFT+U^56–60 

where $UnlI$ is applied to the nl subshell of the Ith atom (e.g., 1s in a single H atom shown in Fig. 2) and σ is the spin index. The occupation matrix, $nnlI,σ$⁠, contains the elements

where, in practice and in this work, each of the v extended states, ψ, for each k point is projected onto atomic orbitals, ϕ, that are obtained during generation of the relevant pseudopotentials.

Like hybrids,^35,76,82,91 the simplified form of DFT+U can only improve FCL errors at the cost of worsening FSL errors.⁷¹ To overcome this limitation, we also incorporate the less widely applied J term^57,92

where the first term in this equation corresponds to Eq. (1) except that U has been replaced by U − J. The J term alone is insufficient to eliminate energetic errors across the full flat plane due to its monotonic dependence on electron count.⁷¹ For s-electron valences, we introduced a modified J (i.e., J′) term

This jmDFT U + J/J′ approach thus has four total coefficients: J and U₁ (U₁ = U − J) values on the n_α + n_β ≤ 1 side and J′ and U₂ (U₂ = U − J′) values on the N = n_α + n_β > 1 side where U is the quantity in Eq. (1). The total expression that applies is

Throughout this work, the FSL belongs to the N = n_α + n_β ≤ 1 side (Fig. 2).

We introduce expressions for nonempirical coefficients and assess their agreement with values obtained from least squares fits of the GGA flat plane errors. In building nonempirical coefficients for model-Hamiltonian inspired terms in jmDFT, coefficients are distinct from those that would be used in DFT+U because the correction to recover piecewise linearity^69,93 differs fundamentally from computed Hubbard U terms.^60,94,95 We note that a single orbital would need to be selected for the modified J′ term as it is written in Eq. (5), which introduces ambiguity into this expression for p and d valence ions. For the p and d valence ions described, we also quantify the elimination of errors for both nonempirical and fitting coefficients in the n_α + n_β ≤ 1 side, whereas we characterize full flat plane errors for the s valence ions. The same ambiguity does not extend, however, to s valence molecules: as long as the ionization process is captured by the occupation definitions, a flat plane may still be identified and corrected as we showed in prior work.^69,71

Spin-polarized DFT calculations were carried out with the plane-wave periodic boundary condition code Quantum-ESPRESSO⁹⁶ v5.1 and used the Perdew-Burke-Ernzerhof (PBE)⁹⁷ GGA as the semilocal DFT functional. Ultrasoft pseudopotentials⁹⁸ (USPPs) were employed for all atoms (H, Li, Be, Na, B, C, N, Cr, and Fe), except in cases (He and Mg) where only PBE norm-conserving pseudopotentials⁹⁹ (NCPPs) were available at the time of the study. Semicore states were included in the valence for select systems (Li, 1s; Na or Mg, 2s/2p; and Cr or Fe, 3s/3p). All pseudopotentials were obtained from the Quantum-ESPRESSO website,¹⁰⁰ and a list is provided in the supplementary material, Table S1. Cutoffs of 80 Ry for the wavefunction and 400 Ry for the charge density were used for the He and Mg NCPPs along with Li, Be, B, and Na. For the remaining USPPs (H, C, N, Cr, and Fe), cutoffs of 30 Ry for the wavefunction and 300 Ry for the density were used, consistent with testing in prior work.⁶⁹ 

To enable comparison of total energies of systems of differing charge, the Martyna-Tuckerman scheme¹⁰¹ was employed in conjunction with a 14.8 Å cubic box. Fractional charge calculations were carried out with manually input band (i.e., orbital) occupations in step sizes of 0.1 electrons across both n_α and n_β, and no symmetry of the wavefunction enforced with electron configurations described in the main text. Total magnetization is not fixed in these calculations, but no smearing is incorporated, leading to final electron configurations matching the input. Deviations from desired spatial spin (i.e., spin contamination), if observed, occur for cases with input n_β > n_α, in which case we report the fully equivalent transposed n_α > n_β case. In all cases, no fewer than 10 and up to 15 unoccupied bands were included. To aid self-consistent field (SCF) convergence, the mixing factor was reduced to 0.4 from its default, and the convergence threshold for the SCF energy error was increased slightly to 9 × 10⁻⁶ Ry.

Exact flat planes were constructed as in prior work⁷¹ with only small differences. For fitting coefficients in this work, PBE ionization energies and electron affinities were used to characterize the flat plane error because the U + J/J′ functional used here has no terms to correct errors at integer electron count. After computing deviations of fractional points from exact flat plane behavior, an in-house MATLAB¹⁰² script was used to obtain coefficients from a nonlinear least squares fitting procedure over the whole flat plane. Spin and ionization states of elements used in flat planes or FCLs are provided in the supplementary material, Table S1. The J′ correction was implemented in a locally modified copy of Quantum-ESPRESSO⁹⁶ v5.1, with the source code available upon request.

Densities (for H and He) were obtained with the pp.x postprocessing tool from PBE and jmDFT calculations in Quantum-ESPRESSO.⁹⁶ Reference densities were obtained with ORCA¹⁰³ v4.0.0.2, using orbital-optimized coupled cluster singles and doubles with perturbative triples [i.e., OOCCSD(T)¹⁰⁴] and the cc-pVTZ¹⁰⁵ basis set. For all calculations, densities were obtained with a 160 × 160 × 160 real space grid resolution in a 14.8 Å box (i.e., distance between points: 0.0925 Å). Radial probability and cumulative distribution functions (RPDF and CDF, respectively) were obtained from these grids with an in-house python script as follows. For the RPDF, the density was summed at increasing distances, r, of grid point locations, $r→k$⁠, to the atom center, $c→$⁠,

where the bin length, Δr, is twice the grid resolution (0.185 Å). The CDF at r is obtained from numerical integration of the RPDF for all i points r_i < r

The jmDFT correction along the FSL by definition contains both the U and J terms of Eq. (3), whereas only the first term is nonzero along the edge FCLs (i.e., FCL+ or FCL0, see Sec. II). Thus, in order to obtain a nonempirical expression of the FSL, our first step is to confirm robust nonempirical estimation^11,69,71,93 of the deviation from linearity^{78–80,86,87,106} along the FCL. For the low order fit in jmDFT with U + J/J′, a constant curvature⁹³ with fractional charge, q, should hold

where the FCL here is between the N-1- and N-electron systems. In the jmDFT functional outlined in Sec. II, Eq. (8) corresponds to U₁ to be applied to the FCL+ in the limit of 100% efficiency expected for atoms.⁶⁹ An analogous expression applies for FCL0 based on the HOMO and LUMO of the N+1- and N-electron system, respectively. This expression may also be obtained by examining a cubic spline expression¹¹ 

where q is the fractional charge that varies between 0 for the N-1-electron system and 1 for the N-electron system, and the N-1-electron LUMO and N-electron HOMO eigenvalues estimate the instantaneous slope. The ΔE term is the negative of the ionization potential (negative of the electron affinity) of the N-electron (N-1-electron) system

This spline was introduced¹¹ and recently confirmed¹⁰⁷ to estimate deviation from linearity without explicit fractional charge calculations. If the maximum deviation holds at q = 0.5, then the contributions from ΔE cancel, and the constant curvature expression⁹³ is recovered. For us to develop nonempirical expressions for J expressions in the context of jmDFT U + J/J′ flat plane recovery, these low order, constant curvature expressions must hold for recovering piecewise linearity along the FCLs.

Equation (9) fails when (i) a higher order fit is necessary, i.e., the point of maximum deviation is at a q value that differs from 0.5, or (ii) the eigenvalues are a poor approximation to the first derivative of the energy at the endpoints. In the first case, our jmDFT approach will also fail if higher order terms contribute significantly to the flat plane error. In cases where this occurs, more flexible definitions of the functional correction outlined in previous work⁷¹ would be necessary. For the second type of errors, the LUMO may be poorly defined, particularly where it corresponds to a zero-electron system (e.g., He²⁺ or H⁺) or systems with zero electrons in the spin manifold (e.g., He⁺ or H). Fractional charge energies for He along the FCL+ and FCL0 are indeed poor fits to the cubic spline in Eq. (9), despite the fact that fractional charge energies appear symmetric around q = 0.5, and a quadratic fit performs well (Fig. 3). Thus, this is a case of the second type of error as both the He²⁺ and He⁺ LUMOs are poor approximations of the first derivatives of the energy (Fig. 3). Derivatives of a fit of fractional charge energies are also slightly unstable near the endpoints due to noise in the fractional charge calculations near the q = 0 endpoint of each FCL (Fig. 3). Similar behavior of the LUMOs and derived curvatures occurs for other light s-block elements, especially H⁺ and Li⁺ (supplementary material, Figs. S1–S5).

To overcome this second type of error, we define a curvature expression with a symmetric form that only employs the HOMO for each FCL and its deviation from ΔE (red dotted lines in Fig. 3). The motivation for this approach can be seen by returning to the cubic spline and noting that the constant curvature approximation is the same as the assumption that the eigenvalue contributions to the expression are equivalent

When one of the eigenvalues is unreliable, this equality will not hold, and the interpolation will be skewed even if the computed energies exhibit quadratic dependence with fractional charge (Fig. 3). Instead, we can replace the LUMO term with a HOMO term for the q → 0 portion of the spline and obtain an expression for the curvature as

This symmetric expression recovers good agreement with the quadratic approximation in both He FCLs and other s-block elements (Fig. 3, additional elements and details provided in the supplementary material, Figs. S1–S5 and Text S1).

Over a large range of s-block (H, Li, Na, He, Be, and Mg), p-block (B, C, N), and d-block (Cr or Fe) elements, good agreement is observed between the constant curvature approximation and fits of the full flat plane data (Fig. 4 and supplementary material, Figs. S1–S11). There is no apparent need to employ the symmetric expression for p-block and d-block elements, including those that approach fully ionized valence electron configurations (e.g., d⁰ Cr⁶⁺ or p⁰ N³⁺, see Fig. 4 and supplementary material, Fig. S11). Conversely, for nearly all s-block cases, some improvement is observed (Fig. 4). Second derivatives of the fractional charge fits are also smoother for the p-block and d-block elements although some second derivatives have a nonzero slope that suggests third order expressions could be beneficial in future jmDFT development (supplementary material, Figs. S6–S11).

The benefit of both the symmetric expression and constant curvature approximations discussed here is that neither requires explicit fractional charge calculation. Therefore, we also suggest a metric based on integer electron properties that can be used to decide when the symmetric expression is needed or constant curvature is sufficient. This metric, m, is the absolute value of the ratio of LUMO and HOMO errors⁶⁹ 

where an analogous expression for this FCL+ metric may be written for other FCLs. The constant curvature approximation holds exactly when m is 1.0, and the symmetric expression is beneficial for m < 1.0 (i.e., all s-block FCLs, supplementary material, Tables S2–S4). Using the m = 1.0 cutoff to select the best expression for the curvature then improves agreement of nonempirical coefficients with fitting values for s-block elements and leaves the p- and d-block values unchanged (Fig. 4 and supplementary material, Fig. S11 and Tables S2–S4). Using the resulting nonempirical curvature values leads to modest errors (0.0–0.04 eV on average) with respect to values obtained from empirical fitting (supplementary material, Tables S2–S4).

Returning to the elemental dependence of curvature, most values span a 5–20 eV range, regardless of angular momentum (Fig. 4 and supplementary material, Fig. S11). As expected,¹⁴ successive ionizations increase the curvature as does increasing nuclear charge for isoelectronic species, i.e., in the same row of the periodic table (Fig. 4 and supplementary material, Fig. S11). The incorporation of semicore states in the valence (e.g., Li, Na, or Mg) has a limited effect on the computed curvature in comparison with other elements with only valence electrons explicitly described (Fig. 4). Overall, sufficiently strong chemical dependence is observed to motivate nonempirical evaluation rather than assuming a universal set of coefficients, but some transferability (e.g., between neighboring oxidation states) should be possible (Fig. 4 and supplementary material, Fig. S11).

We have shown that the nonempirical constant curvature approximation or a modified, symmetric expression reproduce the observed curvature along the FCLs. For our strategy to work, it is also essential that these U values recover piecewise linearity along the FCL. We previously showed for an Fe³⁺ ion that piecewise linearity was recovered exactly at the point anticipated by applying the curvature value obtained in the framework of DFT+U. However, lower correction efficiencies were sometimes observed for transition metal complexes due to hybridization.⁶⁹ Using the nonempirical U values derived in this work, deviations from linearity are eliminated over a set of 10 representative s-block, p-block, and d-block ions (supplementary material, Table S5). Deviations from linearity at q = 0.5 of approximately 1–3 eV for PBE are reduced by over 95% to below 0.05 eV for the FCLs studied (supplementary material, Table S5). The degree of linearity recovered is comparable for cases where the symmetric nonempirical U is the appropriate coefficient and where the constant curvature value is suitable (supplementary material, Table S5).

We have previously demonstrated how to use U or U + J terms to recover piecewise linearity⁶⁹ and the flat plane,⁷¹ respectively. We now introduce an analytical expression for the J coefficient. In the present single site, single subshell case, the U + J correction has the general form

For a single n_α and n_β electron, the expression simplifies further

where n_α and n_β are the individual electron counts and we have neglected any effect of off-diagonal terms on the evaluation of Eq. (14) (i.e., strictly correct for s valence but approximate for higher angular momentum). The FSL is defined as the portion of the flat plane, where n_β = 1− n_α, thus simplifying our expression further to

where the U₁ term is the nonempirical U value obtained from determining the deviation from linearity between the N-1- and N-electron systems [e.g., from Eq. (8)] and J is the analytical coefficient we aim to define here. Energetics along the FSL are affected by both the U and J terms simultaneously, meaning that (i) challenges in U determination for some systems (e.g., H, He⁺) could propagate here and will need to be addressed and (ii) the J term by definition must offset the increased FSL errors caused by the U correction.⁷¹ 

The coefficients we wish to determine are the second derivative of the energy with respect to variation of n_α(n_β) along the FSL (i.e., where n_α + n_β = 1)

where this second derivative should be equal to zero for a functional that exhibits no errors on the FSL. For an N-electron system starting with a single spin-up (i.e., α) electron, this FSL curvature is associated with the instantaneous difference in energetic change for loss of an infinitesimal amount of α electron and gain of an infinitesimal amount of spin down (i.e., β) electron. These terms may be recognized as the spin-up HOMO and spin-down LUMO for the N-electron case, which yields

This expression can also be recognized as a finite difference form of the second derivative along the fractional spin line. Strictly speaking, this expression corresponds to the energetic cost of removing an electron from the α HOMO to infinity and bringing a different electron from infinity to populate the β LUMO. The quantity we are interested in is more directly related to the energetic cost of shifting a fractional electron between spin states without changing the number of electrons. However, we will demonstrate that this difference has no impact on numerical values, unlike curvature coefficients.⁶⁹ 

Since we already have a definition for U₁, we can simplify our expression and define the nonempirical J coefficient as

Before proceeding further to validate this expression, we explore the relationship between the orbital degeneracy and static correlation error (SCE). The SCE of the systems studied in this work can be evaluated as the difference in energy of the system in which an electron is split into half over isoenergetic α and β orbitals vs placed in a single α or β orbital

The SCE corresponds to the extremum along the FSL, which is by definition symmetric in contrast with the curvature expression that need not be (see Sec. IV A). We also define the magnitude of the difference in α HOMO and β LUMO eigenvalues, |Δε|, at the N-electron (n_α = 1, n_β = 0) endpoint

For nearly all density functional approximations (DFAs), this expression will be positive, making the absolute value unnecessary. A comparison of these two quantities in s¹, p¹, and d¹ electron configurations confirms a good correlation between SCE and |Δε| (Fig. 5). The lightest elements have very large eigenvalue deviations and SCE, whereas boron and heavier elements have more modest values for both quantities (Fig. 5). These results are also consistent with prior observations^76,91 that both quantities are relatively small with a semilocal GGA. A single linear correlation between the two quantities holds (R² = 0.99) over ten elements and multiple principle quantum numbers, both with and without semicore states in the valence.

We next directly compare the magnitude of J values obtained from the jmDFT quadratic fit over the full flat plane to nonempirical values for s¹ electron configurations. Fully nonempirical J values are in excellent agreement with fitting values (i.e., within 0.3 eV) for higher Z (i.e., Be⁺, Na, and Mg⁺) elements (Fig. 6 and supplementary material, Table S6). As expected, larger errors are observed for H and He⁺, leading to underestimates of |J| by 4.1 and 4.7 eV, respectively, due to challenges with directly employing LUMO eigenvalues for these two elements (Fig. 6 and see Fig. 3 and supplementary material. Fig. S1). The Li case has an intermediate 1.0 eV error although the relative errors in the nonempirical |J| value for He⁺ and Li are comparable at around 15% and 17%, respectively (Fig. 6 and supplementary material, Table S6).

In Sec. IV A, we showed that our nonempirical J expression depends on both the LUMO of the N-1-electron and N-electron systems, which appeared to be poor approximations to first derivatives for H, He⁺, and Li especially (Fig. 3 and supplementary material, Figs. S1 and S2). We isolate this source of error in the J expression by subtracting a more accurate form of U₁ [e.g., from Eq. (12) or from fitting, U_1,fit] leaving an expression for J that depends only on the HOMO and LUMO of the N-electron system

Selecting the J@U_1,fit significantly improves errors in H (1.5 eV) and Li (0.0 eV) and has no effect on errors for most systems where J was already well predicted (i.e., Na and Mg⁺, Fig. 6 and supplementary material, Table S6). At the same time, somewhat worsened He⁺ and Be⁺ errors suggest that expressions that depend on LUMO eigenvalues will remain problematic in challenging cases.

As for the FCL curvature, we define an expression for J that avoids dependence on any LUMO eigenvalues. Starting from the N-electron configuration with n_α = $12$⁠, n_β = $12$ occupations, we follow an FCL where n_α = n_β toward the N-1-electron system. Along this FCL, the slope of the jmDFT U+J correction is exactly J/2. The slope of the flat-plane error along this FCL is the deviation of the half-filled frontier orbital energy at the n_α = $12$⁠, n_β = $12$ configuration from the negative of the ionization potential of the N-electron system, ΔE,

This expression produces by far the best agreement with the fitting coefficients for challenging H and He⁺ cases and performs comparably to other nonempirical expressions for the remaining elements but requires a fractional spin calculation (Fig. 6).

The preferred expression for determining a nonempirical J coefficient will depend on the use case. The most robust expression [Eq. (23)] requires a fractional spin calculation, which is not widely available with all electronic structure codes and methods. Where fractional spin calculations are to be avoided, an alternative approach is to use the other expressions for J to obtain an interpolation of errors along the fractional spin line

only from properties obtained from the N-electron system, in analogy to the expressions employed for FCL interpolation.¹¹ For cases where the nonempirical U and J expressions are in good agreement with the fit, interpolated values are in excellent agreement with the explicit fractional spin calculations (supplementary material, Fig. S12).

Using these nonempirical coefficients, it is useful to next determine if the static correlation error along the FSL can be minimized and if degeneracy of the frontier orbitals can be restored. We incorporate the nonempirical J term obtained using the symmetric form of U (i.e., J@U_1,symm) alongside the symmetric, nonempirical U₁ term. U₁ alone would normally increase SCE, but the J term should mitigate this error, recovering invariant energies and frontier MO eigenvalues across the FSL. Indeed, we observe this behavior for a representative s¹ system, Mg⁺, with semilocal DFT SCE of over 0.3 eV reduced to less than 0.05 eV with nonempirical jmDFT (Fig. 7). At the same time, the frontier eigenvalues become near degenerate across the range of n_α–n_β values with jmDFT in contrast to semilocal DFT where they differ in energy by as much as 1.5 eV (Fig. 7). These observations hold in Be⁺ and Li as well, reducing SCE as large as 0.6 eV by nearly an order of magnitude and recovering equivalence of the α and β frontier orbitals (supplementary material, Figs. S13 and S14). The use of DFT+U alone with the same U₁ values that eliminate FCL errors instead increases semilocal DFT errors by an order of magnitude in both total and orbital energies (supplementary material, Figs. S13 and S14).

Over all s¹ systems, we confirm that nonempirical J coefficients selected in conjunction with the nonempirical, symmetric U₁ (i.e., J@U_1,symm) dramatically reduce differences in frontier orbital energies (i.e., |Δε|) (Fig. 8). The differences for Be⁺, Na⁺, and Mg⁺ are essentially eliminated, whereas the cases that had larger sensitivities to LUMO energetics (i.e., H and He⁺) still have a significant residual |Δε|. Nearly 90% of this error is eliminated across all systems under study, including H or He⁺ (Fig. 8). We note that the error metric chosen here depends on the LUMO, which further complicates assessment. Self-consistent adjustment of J until frontier MO eigenvalues become degenerate could be an alternative to the nonempirical approach we have pursued.

The nonempirical J expression also holds for higher angular momentum electrons (i.e., p¹: B, C⁺, and N²⁺ and d¹: Cr⁵⁺). We first obtained the value of U₁ + J from fitting to the PBE FSL error. The resulting value is in very good agreement with the nonempirical U₁ + J obtained from eigenvalues (supplementary material, Fig. S15 and Table S7). After subtracting a nonempirical U₁ value to obtain the fitting J value, differences between the fit and nonempirical |J| coefficients for the p-block and d-block ions are within 0.1 eV (Fig. 9). These nonempirical coefficients are better predictions of the fitting coefficient in p- and d-electron systems than in the s-electron ions, consistent with observations (see Sec. IV A) on U evaluation (see Fig. 6).

We next considered the extent to which SCE elimination is achieved at the nonempirical coefficient values. Because the nonempirical and fitting coefficients are in such good agreement, the SCE elimination reflects the efficiency of the jmDFT U + J correction on these systems. In all four p¹ and d¹ cases studied here, the semilocal DFT SCE is quite modest (i.e., 0.2–0.4 eV) but increases dramatically (i.e., 2.5–4.5 eV) when the nonempirical U coefficient is applied (supplementary material, Table S7). In representative p-electron B and d-electron Cr⁵⁺, the nonempirical U + J correction eliminates most of the SCE introduced by the U-only correction, but typically has residual SCE slightly above the PBE value (0.3–0.7 eV, Fig. 9 and supplementary material, Table S7). Over p-electron (B, C⁺, N²⁺) and d-electron (Cr⁵⁺) elements, we are still able to eliminate SCE fully by increasing the applied |J| value between 1 and 4 eV above its nonempirical value (Fig. 9 and supplementary material, Fig. S16 and Table S7).

The discrepancy between the J value that eliminates SCE and the nonempirical value arises from two sources. First, 100% efficiency in jmDFT SCE elimination should occur at 0.25 eV/eV of |J| applied because the correction along the FSL is of the order Jn_αn_β, n_α = n_β = 0.5 at the point where SCE is evaluated, and thus the correction applied is of the order J/4. However, a 0.25 eV/eV |J| slope is observed only for B, with lower values especially for Cr⁵⁺ of 0.18 eV/eV |J| (Fig. 9 and supplementary material, Table S7). Second, applying only DFT + U breaks the symmetry of the orbitals, causing the Tr[n_αn_β] term to be zero, but larger |J| values (i.e., |J| > 1–2 eV) return the symmetry of the orbitals (supplementary material, Fig. S17). The change in orbital occupancy leads to an effectively higher SCE state (i.e., if its energy were extrapolated back to |J| = 0 eV) but one in which the U + J functional then eliminates SCE (supplementary material, Fig. S16). This observation is consistent with the fact that the U-only SCE is below the y-intercept of the SCE vs |J| line for higher |J| values by 0.3 (for B) to 0.8 (for N²⁺) eV (Fig. 9 and supplementary material, Fig. S16). A practical strategy to overcome this limitation is to obtain the SCE or |Δε| (i.e., to avoid fractional spin calculations) at several values of |J| and use the best-fit line to identify the |J| value that eliminates the error (Fig. 9 and supplementary material, Table S7 and Fig. S16). This strategy would work for all p- and d-electron systems studied here because the SCE vs |J| relationship is strictly linear for all |J| above around 2–3 eV (Fig. 9 and supplementary material, Fig. S16).

Although the FSL is on the n_α + n_β ≤ 1 side of the flat plane, we can obtain nonempirical U₂ and J′ coefficients for n_α + n_β ≥ 1 by asymptotically approaching the FSL from the other direction. The expression for these n_α + n_β ≥ 1 coefficients is then

where U₂ was already defined as

and J′ is thus

When fitting the full flat plane, it is not necessary to enforce U₂ + J′ = U₁ + J, but eliminating small differences would eliminate unphysical discontinuities near the FSL (supplementary material, Table S8). Because the J′ expression only depends on occupied orbitals, we expect good agreement between the fit and nonempirical coefficients. Indeed, all nonempirical J′ values disagree by no more than 0.4 eV with the fit values [mean absolute error (MAE): 0.18 eV] across s-, p-, and d-electron systems (supplementary material, Table S8). Errors are comparable for all elements and angular momenta, and they are a small percentage of the J′ value (∼1%–3%, see supplementary material, Table S8). The evaluation of J′ solely from occupied orbitals thus reinforces the predictive capabilities of the nonempirical expressions.

We have derived a full set of coefficients suitable for the jmDFT U + J/J′ correction, i.e., U₁ in Eq. (8), J in Eq. (19), U₂ in Eq. (26), and J′ in Eq. (27), with a symmetric U₁ and J obtained from Eqs. (12) and (22) for s-electron cases, as described in Sec. IV A (supplementary material, Tables S2–S9). We now evaluate whether these nonempirical coefficients are effective at eliminating the flat plane error across a range of angular momentum values. The mean absolute errors (MAEs) for the PBE N-1- to N-electron portion of the flat plane are at least 0.5 eV for most systems but smallest for Li and Na (∼0.3 eV) and largest for Cr⁵⁺ and He⁺ at around 1.5 eV (Fig. 10). Increasing ionization generally increases these flat plane errors, whereas increasing atomic number has a more limited effect except for the distinction between H/He⁺ and heavier ions (Fig. 10). Across most systems, nonempirical coefficients are in good agreement with values obtained from low order fits, suggesting that if the latter performs well to eliminate errors,⁷¹ the former should as well (supplementary material, Tables S8–S9). This observation is not limited to the ions that are the focus of this work but also holds true in molecules for which we had only previously applied jmDFT with empirical coefficients⁷¹ (supplementary material, Fig. S18).

Indeed, both sets of coefficients in jmDFT lead to uniformly reduced flat plane MAEs by at least 75% to as low as around 0.05 eV for the best performing s-electron systems (Li, Na, Mg, see Fig. 10 and supplementary material, Table S9). Here, we apply the symmetric form of the nonempirical U₁ and J for the s-electron systems to avoid using the LUMO eigenvalue expressions in cases where it performed unreliably. In H and He⁺, for which coefficient estimation was identified as challenging (see Sec. IV A), fitting results improve slightly over the nonempirical values (Fig. 10). At the same time, the large PBE errors for H/He⁺ mean that the percent reduction in error by jmDFT is actually comparable to that for the other s-electron ions (Fig. 10). Comparing p¹ and d¹ ions, we observe that large PBE errors are reduced significantly in all systems equivalently with nonempirical and fitting coefficients (Fig. 10). We had previously noted (see Sec. IV B) that differences in the PBE or DFT+U ground state led to systematic, numerical underestimation of the needed J values to eliminate FSL errors, whereas the FCL errors were straightforward (see Sec. IV A) to eliminate. Indeed, more residual errors remain on the FSL than on the FCL for these systems (supplementary material, Table S10).

In addition to errors on the n_α + n_β ≤ 1 side of the plane, determination of full flat-plane errors for s-electron systems is also possible (supplementary material, Table S9). In Mg⁺, good recovery of the full flat plane is observed on both sides using the U + J/J′ correction (Fig. 10). Although errors are overall reduced with respect to PBE, jmDFT performs better for cases where the N- to N+1-electron side corresponds to a neutral species (e.g., He⁺) and less well for the cases that correspond to an anion (e.g., H⁻, supplementary material, Table S9). Examining occupations as fractional electrons are added to H reveals that not all fractional electrons are added to the appropriate s orbital within PBE, and the jmDFT correction worsens the discrepancy between expected and observed occupations (supplementary material, Fig. S19). As outlined in Sec. II, when occupations are no longer a good indicator of electron count, this form of jmDFT corrections is expected to fail, and a different approach is needed.

Given the renewed focus^108–111 on the relationship between the density and energetic error,^9,43,112 especially relevant for SIE-prone systems (e.g., transition metals^{59,74–76,113,114}), we revisit the effect of a jmDFT U + J/J′ correction on densities. We compare representative PBE and jmDFT H⁻ and He densities to more accurate wavefunction theory densities [i.e., OOCCSD(T), see Sec. III]. We compared the radial probability and cumulative distribution functions (i.e., RPDFs and CDFs), g(r) and G(r), which quantify the instantaneous and cumulative probability, respectively, of finding electron density at a specific distance from the nucleus. We focus on the H⁻ anion, which is representative of other anion N+1-electron systems (i.e., Li⁻, Na⁻, B⁻) that have positive HOMO eigenvalues consistent with weakly bound electrons in semilocal DFT.^{14,42,115,116} The PBE CDF evaluated 10 bohr from the nucleus for H⁻ confirms that PBE only partially binds the 1s² electrons (1.83e) in contrast with the expected value of 2e and the CCSD(T) result (1.99e). Applying a nonempirical U + J/J′ jmDFT correction worsens the discrepancy, with a CDF of 1.72e evaluated at 10 bohrs. Comparing the RPDFs, we observe CCSD(T) density to have a single peak at 2.0 bohrs, whereas both PBE and jmDFT have both a slightly more localized peak at 1.5 bohrs and a second, fully delocalized 15 bohr peak for the partially unbound electron density (Fig. 11). This second distant peak that is already present in PBE is made slightly taller by jmDFT, consistent with observations from the CDFs (Fig. 11). This observation is not surprising as the jmDFT potential can be expected to destabilize the positive H⁻ anion HOMO further. It is also consistent with our prior observations that the underestimated electron affinity of H⁻ from PBE was worsened slightly by jmDFT.⁷¹ 

We also compare densities obtained for the isoelectronic, neutral He N+1-electron system, which instead has a negative HOMO eigenvalue that is representative of the neutral or cationic N+1-electron ions we studied. The CDF evaluated at 10 bohrs for He with PBE (1.99e) and jmDFT (1.99e) are in good agreement with the CCSD(T) reference value (1.99e). All three methods also exhibit comparable RPDFs with a single peak at 1 bohr for He (Fig. 11). Thus, we may conclude that jmDFT will not worsen good PBE densities, but the potential form we have proposed here is unlikely to improve cases where semilocal DFT densities are poor. To overcome this shortcoming, we would need to develop a multipronged correction approach for challenging cases such as H⁻ that applies a potential to localize the density first and then applies energetic corrections. Such approaches are currently under development in our lab.

We have developed a nonempirical approach to eliminating deviations from exact behavior in semilocal DFT. We took a model-Hamiltonian-inspired approach, wherein physically grounded but no-cost algebraic expressions are incorporated into semilocal DFT functionals to correct errors without adding computational cost. We developed nonempirical expressions for each coefficient in the jmDFT approach with the U + J/J′ correction. The U term was first confirmed to be determined by expressions used in prior work that approximate curvature along the fractional charge line. However, we noted potential challenges in employing expressions that depended on unoccupied molecular orbitals for specific cases (e.g., H and He⁺). We introduced both an alternate expression that depended only on occupied eigenvalues and a metric derived from integer-electron properties to determine when the alternate expression should be employed. This approach provided better estimates when compared with explicit fractional charge results. We introduced estimates of the J coefficients that are related to the fractional spin error. We showed that these quantities could be obtained fully from integer charge eigenvalues. Such expressions could in turn be used to estimate the fractional spin error in cases where explicit fractional occupancy calculations were infeasible.

We demonstrated elimination of both fractional charge and fractional spin errors at nonempirical coefficient values. In the case of applying J to eliminate the fractional spin error, we noted that changes in the occupation of orbitals meant that PBE or DFT+U eigenvalues would underestimate the needed |J| but that systematically increasing |J| values eliminated the fractional spin error in s¹, p¹, and d¹ systems. We demonstrated that nonempirical coefficients eliminated the PBE flat plane error, as judged through reduction in energetic MAEs, as well as the coefficients obtained from fitting the original PBE error. In all cases, errors on the FCL were reduced significantly (∼95%) without the increase in FSL errors that would occur for a standard approach (e.g., DFT+U or hybrids). Further improvement in most cases could be achieved through an iterative approach, e.g., as shown for the |J| adjustments, to account for evolving electronic structure when both U and J terms are incorporated. Finally, we noted that outstanding challenges in eliminating the flat plane error remain when electrons do not localize well to the targeted subshell, especially in weakly bound anions. We observed in these cases that jmDFT worsened PBE densities and occupations. In nonpathological cases, jmDFT and PBE had comparably accurate densities, as judged against WFT references. Alternative approaches that first stabilize the anion prior to energetic correction will be a necessary step forward to improve both density and energetic properties. Further development of flexible jmDFT expressions will also be necessary to enable the practical application of this approach to transition metal complexes and other larger molecules. Such efforts are currently underway in our lab.

See supplementary material for the list of pseudopotentials used; deviations from piecewise linearity, slopes, and curvature of H, Li, Be, Na, Mg, B, C, N, Cr, and Fe; details on curvature for p- and d-electron valence ions; alternative nonempirical expressions for U₁ and U₂; symmetric and constant curvature spline result comparison to fits for s, p, and d curvatures; nonempirical fractional charge line recovery results; nonempirical estimation of J values; FSL interpolation estimates; energetics and frontier MO eigenvalues on FSL with corrections for Be⁺ and Li; details of U + J from fitting and nonempirical results for FSL recovery on p and d systems; FSL error for C⁺ and N²⁺ with J varied; comparison of occupations and FSL errors for B and Cr⁵⁺ with J varied; comparison of J/J′ values obtained from different expressions; results on the H₂ molecule; flat plane error estimations; and all total energies and fitting scripts.

The authors acknowledge support by the Department of Energy under Grant No. DE-SC0018096. H.J.K. holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund. The authors thank Adam H. Steeves for providing a critical reading of the manuscript.