Density functional theory (DFT) is widely applied to both molecules and materials, but well known energetic delocalization and static correlation errors in practical exchange-correlation approximations limit quantitative accuracy. Common methods that correct energetic delocalization errors, such as the Hubbard U correction in DFT+U or Hartree-Fock exchange in global hybrids, do so at the cost of worsening static correlation errors. We recently introduced an alternate approach [Bajaj et al., J. Chem. Phys. 147, 191101 (2017)] known as judiciously modified DFT (jmDFT), wherein the deviation from exact behavior of semilocal functionals over both fractional spin and charge, i.e., the so-called flat plane, was used to motivate functional forms of second order analytic corrections. In this work, we introduce fully nonempirical expressions for all four coefficients in a DFT+U+J-inspired form of jmDFT, where all coefficients are obtained only from energies and eigenvalues of the integer-electron systems. We show good agreement for U and J coefficients obtained nonempirically as compared with the results of numerical fitting in a jmDFT U+J/J′ correction. Incorporating the fully nonempirical jmDFT correction reduces and even eliminates the fractional spin error at the same time as eliminating the energetic delocalization error. We show that this approach extends beyond s-electron systems to higher angular momentum cases including p- and d-electrons. Finally, we diagnose some shortcomings of the current jmDFT approach that limit its ability to improve upon DFT results for cases such as weakly bound anions due to poor underlying semilocal functional behavior.

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 error12–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,19 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 pursued29–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 error36–41 or other schemes to localize orbitals.9,42,43 Generalizations44 to Kohn-Sham DFT, e.g., global or range-separated hybrids45–55 and DFT+U,56–60 fall between these two limits. The efficacy of these approaches in improving materials and molecular properties has been explained4,7,11,20,59–70 by noting that they generally reduce DE. A number of variants71–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 density74,75 or static correlation error (SCE),76 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) approach71 to simultaneously correct static correlation and delocalization errors. A limitation of the initial jmDFT effort71 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 introduced71 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,71 we employ the flat plane78–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 charge56–59 and unchanged energy81,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 methods83,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 discontinuity79,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).

FIG. 1.

A typical flat plane constructed from interpolated energies (in electron volts) as a function of (nα, nβ). The exact energies (black circles) are (i) the N-1-electron system, E(0, 0); (ii) the N-electron system, E(1, 0); and (iii) the N+1-electron system, E(1, 1). The FSL is highlighted in blue and connects isoenergetic E(1, 0) and E(0, 1) configurations. The FCL+ represents the ionization potential of the N-electron system, the FCL0 represents the electron affinity of the N-electron system, and both are shown in dark red.

FIG. 1.

A typical flat plane constructed from interpolated energies (in electron volts) as a function of (nα, nβ). The exact energies (black circles) are (i) the N-1-electron system, E(0, 0); (ii) the N-electron system, E(1, 0); and (iii) the N+1-electron system, E(1, 1). The FSL is highlighted in blue and connects isoenergetic E(1, 0) and E(0, 1) configurations. The FCL+ represents the ionization potential of the N-electron system, the FCL0 represents the electron affinity of the N-electron system, and both are shown in dark red.

Close modal

Analysis of experimental ionization potentials of gas phase ions88 revealed89 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,71 we focus on electron configurations in which an initially empty valence shell in the N-1-electron ion (e.g., singlet s0) is populated first by an nα electron to form the N-electron ion (e.g., doublet s1) followed by an nβ electron to form the N+1-electron ion (e.g., singlet s2). 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 detail71,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 p0 or d0, N-electron doublet p1 or d1, and N+1-electron singlet p2 or d2. The singlet N+1-electron configurations are no longer ground states as they were in the s valence cases but are well defined90 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 N2+ for p1 ions and Cr5+ for d1 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 suggested71 that low order corrections with few coefficients could be used to correct the majority of the flat plane error. These observations led us to compare71 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.

FIG. 2.

Deviations of the PBE GGA from exact flat plane behavior (Edev, in electron volts) vs electron count (nα, nβ) after alignment of all endpoints in the flat plane to the PBE results. The representative case of H is shown, where (0, 0) is H+, (1, 0) is H atom, and (1, 1) is H. The FSL is highlighted in dark blue, and the FCL+ and FCL0 are both shown in dark red. Values are also projected onto the surface, as indicated by the inset color bar.

FIG. 2.

Deviations of the PBE GGA from exact flat plane behavior (Edev, in electron volts) vs electron count (nα, nβ) after alignment of all endpoints in the flat plane to the PBE results. The representative case of H is shown, where (0, 0) is H+, (1, 0) is H atom, and (1, 1) is H. The FSL is highlighted in dark blue, and the FCL+ and FCL0 are both shown in dark red. Values are also projected onto the surface, as indicated by the inset color bar.

Close modal

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.71 These forms have close connection to the widely applied, rotationally invariant and simplified form of DFT+U56–60 

EDFT+U=EDFT+12I,σnlUnlITrnnlI,σ1nnlI,σ,
(1)

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

nmmI,σ=k,vψk,vσϕmIϕmIψk,vσ,
(2)

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.71 To overcome this limitation, we also incorporate the less widely applied J term57,92

EDFT+U+J=EDFT+12I,σnl(UnlIJnlI)TrnnlI,σ1nnlI,σ+12I,σnlJnlITrnnlI,σnnlI,σ,
(3)

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

EJ=12I,σnlJ  nlITrnnlI,σ1nnlI,σ1.
(4)

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

EU+J/J=12I,σnl(UnlIJnlI)TrnnlI,σ1nnlI,σ+12I,σnlJnlITrnnlI,σnnlI,σifN112I,σnl(UnlIJ  nlI)TrnnlI,σ1nnlI,σ+12I,σnlJ  nlITrnnlI,σ1nnlI,σ1ifN>1.
(5)

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 linearity69,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-ESPRESSO96 v5.1 and used the Perdew-Burke-Ernzerhof (PBE)97 GGA as the semilocal DFT functional. Ultrasoft pseudopotentials98 (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 pseudopotentials99 (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,100 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.69 

To enable comparison of total energies of systems of differing charge, the Martyna-Tuckerman scheme101 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−6 Ry.

Exact flat planes were constructed as in prior work71 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 MATLAB102 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-ESPRESSO96 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.96 Reference densities were obtained with ORCA103 v4.0.0.2, using orbital-optimized coupled cluster singles and doubles with perturbative triples [i.e., OOCCSD(T)104] and the cc-pVTZ105 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, rk, to the atom center, c,

g(r)=k,r12Δr|crk|<r+12Δrρ(rk)ΔVΔr,
(6)

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 ri < r

G(r)=0rg(s)ds=ri<rg(ri)Δr.
(7)

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 estimation11,69,71,93 of the deviation from linearity78–80,86,87,106 along the FCL. For the low order fit in jmDFT with U + J/J′, a constant curvature93 with fractional charge, q, should hold

U1=2Eq2=εNHOMOεN1LUMO,
(8)

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 U1 to be applied to the FCL+ in the limit of 100% efficiency expected for atoms.69 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 expression11 

Edev(q)=E(q)ΔEq=[(εN1LUMOΔE)(1q)+(ΔEεNHOMO)q]q(1q),
(9)

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

ΔE=E(N)E(N-1).
(10)

This spline was introduced11 and recently confirmed107 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 expression93 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 work71 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., He2+ 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 He2+ 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).

FIG. 3.

(top) Deviation from piecewise linearity (Edev, in electron volts) for He from fractional charge calculations along the FCL+ (He2+ to He+, green circles) and FCL0 (He+ to He, blue circles) with quadratic fit dashed line shown in the same color. A cubic spline (orange dotted line) and symmetrized form (red dotted line) are shown along with the first derivative from the endpoint eigenvalues (pink arrows). The midpoint is shown in a dashed black vertical line for both plots. (bottom) Second derivatives of Edev quadratic fit, cubic and symmetric splines taken at the midpoint are all shown as horizontal lines in the same style as the top pane. The second derivative of a fifth order fit of the Edev at top is shown for the 0.2–0.8e range in circle symbols.

FIG. 3.

(top) Deviation from piecewise linearity (Edev, in electron volts) for He from fractional charge calculations along the FCL+ (He2+ to He+, green circles) and FCL0 (He+ to He, blue circles) with quadratic fit dashed line shown in the same color. A cubic spline (orange dotted line) and symmetrized form (red dotted line) are shown along with the first derivative from the endpoint eigenvalues (pink arrows). The midpoint is shown in a dashed black vertical line for both plots. (bottom) Second derivatives of Edev quadratic fit, cubic and symmetric splines taken at the midpoint are all shown as horizontal lines in the same style as the top pane. The second derivative of a fifth order fit of the Edev at top is shown for the 0.2–0.8e range in circle symbols.

Close modal

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

εN1LUMOΔE=ΔEεNHOMO.
(11)

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

U1,symm=2Eq22εNHOMOΔE.
(12)

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., d0 Cr6+ or p0 N3+, 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).

FIG. 4.

Curvature values (i.e., |U|, in electron volts) obtained from a quadratic least squares fit of fractional charge calculations (gray bars) compared to symmetrized U1 values from Eq. (12) (symm., red horizontal lines) and the constant curvature U1 values from Eq. (8) (CC, blue horizontal lines). Results are shown as two panes with FCL+ (top) and FCL0 (bottom) values for s-electron (H, Li, Na, He+, Be+, and Mg+), p-electron (B, C+, and N2+), and d-electron (Cr5+) ionizations, where the x-axis label corresponds to the N-electron species.

FIG. 4.

Curvature values (i.e., |U|, in electron volts) obtained from a quadratic least squares fit of fractional charge calculations (gray bars) compared to symmetrized U1 values from Eq. (12) (symm., red horizontal lines) and the constant curvature U1 values from Eq. (8) (CC, blue horizontal lines). Results are shown as two panes with FCL+ (top) and FCL0 (bottom) values for s-electron (H, Li, Na, He+, Be+, and Mg+), p-electron (B, C+, and N2+), and d-electron (Cr5+) ionizations, where the x-axis label corresponds to the N-electron species.

Close modal

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 errors69 

m=absεN1LUMOΔEΔEεNHOMO,
(13)

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,14 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 Fe3+ 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.69 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 linearity69 and the flat plane,71 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

EU+J=12σU1  Trn  σ1n  σ+12σJ    Trn  σn  σ.
(14)

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

EU+J=12U1  nα(1nα)+(nβ(1nβ))+J(nαnβ),
(15)

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

EU+J=12U1  nα(1nα)+((1nα)(nα))+J(nα(1nα))=(U1  +J)(nα(1nα)),
(16)

where the U1 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.71 

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)

2Enα2=(U1+J),
(17)

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

U1  +JεNHOMOεNLUMO.
(18)

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

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

JεNHOMOεNLUMOεNHOMOεN1LUMO=εN1LUMOεNLUMO.
(19)

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

SCE=E(nα=12,nβ=12)E(nα=1,nβ=0).
(20)

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

Δε=absεNHOMOεNLUMO.
(21)

For nearly all density functional approximations (DFAs), this expression will be positive, making the absolute value unnecessary. A comparison of these two quantities in s1, p1, and d1 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 observations76,91 that both quantities are relatively small with a semilocal GGA. A single linear correlation between the two quantities holds (R2 = 0.99) over ten elements and multiple principle quantum numbers, both with and without semicore states in the valence.

FIG. 5.

Static correlation error (SCE, in electron volts), as judged through the difference in energy between a single electron in a spin up or down orbital vs half the electron in spin up and half in spin down orbitals, vs |Δε| in eV, the unsigned difference in frontier α and β orbitals for the spin up-only configuration. Representative s block (blue circles), p block (red triangles), and d block (green squares) elements are shown, as annotated in the inset, and all points fall on a single line with R2 = 0.997.

FIG. 5.

Static correlation error (SCE, in electron volts), as judged through the difference in energy between a single electron in a spin up or down orbital vs half the electron in spin up and half in spin down orbitals, vs |Δε| in eV, the unsigned difference in frontier α and β orbitals for the spin up-only configuration. Representative s block (blue circles), p block (red triangles), and d block (green squares) elements are shown, as annotated in the inset, and all points fall on a single line with R2 = 0.997.

Close modal

We next directly compare the magnitude of J values obtained from the jmDFT quadratic fit over the full flat plane to nonempirical values for s1 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).

FIG. 6.

Values of |J| in electron volts obtained from a second-order fit of the full flat plane (gray bars), from only the N-electron and N-1-electron LUMO eigenvalues [i.e., J(1, 0) shown in red bars] or using the U obtained from a quadratic fit [i.e., J(1, 0)@Ufit, shown in green bars], and an expression for J obtained from the occupied frontier MO eigenvalues of the (1/2, 1/2) system [i.e., J(1/2, 1/2), blue bars].

FIG. 6.

Values of |J| in electron volts obtained from a second-order fit of the full flat plane (gray bars), from only the N-electron and N-1-electron LUMO eigenvalues [i.e., J(1, 0) shown in red bars] or using the U obtained from a quadratic fit [i.e., J(1, 0)@Ufit, shown in green bars], and an expression for J obtained from the occupied frontier MO eigenvalues of the (1/2, 1/2) system [i.e., J(1/2, 1/2), blue bars].

Close modal

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 U1 [e.g., from Eq. (12) or from fitting, U1,fit] leaving an expression for J that depends only on the HOMO and LUMO of the N-electron system

J@U1,fitεNHOMOεNLUMOU1,fit.
(22)

Selecting the J@U1,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,

J=2ΔEεNHOMOnα=12,nβ=12.
(23)

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

Edev=εNLUMOεNHOMOnα(1nα),
(24)

only from properties obtained from the N-electron system, in analogy to the expressions employed for FCL interpolation.11 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@U1,symm) alongside the symmetric, nonempirical U1 term. U1 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 s1 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 U1 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).

FIG. 7.

Mg+ energetics along the FSL (top, green circles) and frontier α (black open circles) and β (red open circles) molecular orbital energetics (bottom), both in electron volts . The left panes are the results for PBE, and the right panes show the result for jmDFT with U1 = 6.86 eV and J = −8.30 eV.

FIG. 7.

Mg+ energetics along the FSL (top, green circles) and frontier α (black open circles) and β (red open circles) molecular orbital energetics (bottom), both in electron volts . The left panes are the results for PBE, and the right panes show the result for jmDFT with U1 = 6.86 eV and J = −8.30 eV.

Close modal

Over all s1 systems, we confirm that nonempirical J coefficients selected in conjunction with the nonempirical, symmetric U1 (i.e., J@U1,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.

FIG. 8.

The unsigned difference in frontier α and β orbitals for a spin up-only configuration, |Δε|, in electron volts for s1 atoms and ions. The PBE result is shown in red bars, and the jmDFT result with nonempirical U1 and J coefficients (n.e. U + J) is shown in green bars, as indicated in the inset legend.

FIG. 8.

The unsigned difference in frontier α and β orbitals for a spin up-only configuration, |Δε|, in electron volts for s1 atoms and ions. The PBE result is shown in red bars, and the jmDFT result with nonempirical U1 and J coefficients (n.e. U + J) is shown in green bars, as indicated in the inset legend.

Close modal

The nonempirical J expression also holds for higher angular momentum electrons (i.e., p1: B, C+, and N2+ and d1: Cr5+). We first obtained the value of U1 + J from fitting to the PBE FSL error. The resulting value is in very good agreement with the nonempirical U1 + J obtained from eigenvalues (supplementary material, Fig. S15 and Table S7). After subtracting a nonempirical U1 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).

FIG. 9.

Values of |J| (in electron volts ) obtained from a fit to the FSL (fit, red bars) or from the eigenvalue expression [J(1, 0), red bars] for 2p1 B, C+, and N2+ and 3d1 Cr5+ (left pane). Elimination of the static correlation error (SCE, in electron volts ) as judged through deviation of energies at the midway point along the FSL for 2p1 B (middle pane) and 3d1 Cr5+ (right pane) with fixed nonempirical U applied and J varied. Varied J values (gray circles) are shown with a best fit line (gray dashes) along with the original PBE value (horizontal dashed green line), the nonempirical J value (blue square), and the U-only result (red circle). The x-axis of J values is inverted with more negative numbers shown at right.

FIG. 9.

Values of |J| (in electron volts ) obtained from a fit to the FSL (fit, red bars) or from the eigenvalue expression [J(1, 0), red bars] for 2p1 B, C+, and N2+ and 3d1 Cr5+ (left pane). Elimination of the static correlation error (SCE, in electron volts ) as judged through deviation of energies at the midway point along the FSL for 2p1 B (middle pane) and 3d1 Cr5+ (right pane) with fixed nonempirical U applied and J varied. Varied J values (gray circles) are shown with a best fit line (gray dashes) along with the original PBE value (horizontal dashed green line), the nonempirical J value (blue square), and the U-only result (red circle). The x-axis of J values is inverted with more negative numbers shown at right.

Close modal

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 p1 and d1 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 Cr5+, 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+, N2+) and d-electron (Cr5+) 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 Cr5+ 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 N2+) 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 U2 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

U2  +JεNHOMOεNLUMO,
(25)

where U2 was already defined as

U2=εN+1HOMOεNLUMO,
(26)

and J′ is thus

J=εNHOMOεN+1HOMO.
(27)

When fitting the full flat plane, it is not necessary to enforce U2 + J′ = U1 + 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., U1 in Eq. (8), J in Eq. (19), U2 in Eq. (26), and J′ in Eq. (27), with a symmetric U1 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 Cr5+ 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,71 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 coefficients71 (supplementary material, Fig. S18).

FIG. 10.

(a) Example of nonempirical jmDFT (U + J/J′) recovery of the Mg+ flat plane with energies in electron volts and the Mg2+ relative energy set to zero. The x- and y-axis values correspond to number of spin up (nα) and spin down (nβ) 3s electrons, respectively, and charge states of Mg2+ are annotated on the plot. (b) The mean absolute error (MAE, in electron volts) over the N-1- to N-electron portion of the flat plane for ten ions from PBE (red bars), nonempirical jmDFT (U + J) coefficients (n.e. U + J, blue bars), and the result from fit U + J coefficients (green horizontal lines).

FIG. 10.

(a) Example of nonempirical jmDFT (U + J/J′) recovery of the Mg+ flat plane with energies in electron volts and the Mg2+ relative energy set to zero. The x- and y-axis values correspond to number of spin up (nα) and spin down (nβ) 3s electrons, respectively, and charge states of Mg2+ are annotated on the plot. (b) The mean absolute error (MAE, in electron volts) over the N-1- to N-electron portion of the flat plane for ten ions from PBE (red bars), nonempirical jmDFT (U + J) coefficients (n.e. U + J, blue bars), and the result from fit U + J coefficients (green horizontal lines).

Close modal

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 U1 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 p1 and d1 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 focus108–111 on the relationship between the density and energetic error,9,43,112 especially relevant for SIE-prone systems (e.g., transition metals59,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 1s2 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.71 

FIG. 11.

RPDF density profile, γ(r) (e/bohr), for He (top) and He (bottom) obtained with CCSD(T)/cc-pVTZ (gray squares and line), PBE (blue circles and line), and jmDFT (green triangles and line). The van der Waals radii are shown as vertical dashed black lines.

FIG. 11.

RPDF density profile, γ(r) (e/bohr), for He (top) and He (bottom) obtained with CCSD(T)/cc-pVTZ (gray squares and line), PBE (blue circles and line), and jmDFT (green triangles and line). The van der Waals radii are shown as vertical dashed black lines.

Close modal

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 s1, p1, and d1 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 U1 and U2; 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 N2+ with J varied; comparison of occupations and FSL errors for B and Cr5+ with J varied; comparison of J/J′ values obtained from different expressions; results on the H2 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.

1.
A. D.
Becke
,
J. Chem. Phys.
140
,
18A301
(
2014
).
2.
K.
Burke
,
J. Chem. Phys.
136
,
150901
(
2012
).
3.
H. S.
Yu
,
S. L.
Li
, and
D. G.
Truhlar
,
J. Chem. Phys.
145
,
130901
(
2016
).
4.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
J. Chem. Phys.
125
,
201102
(
2006
).
5.
A.
Ruzsinszky
,
J. P.
Perdew
,
G. I.
Csonka
,
O. A.
Vydrov
, and
G. E.
Scuseria
,
J. Chem. Phys.
126
,
104102
(
2007
).
6.
R.
Haunschild
,
T. M.
Henderson
,
C. A.
Jiménez-Hoyos
, and
G. E.
Scuseria
,
J. Chem. Phys.
133
,
134116
(
2010
).
7.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Science
321
,
792
(
2008
).
8.
T.
Schmidt
and
S.
Kümmel
,
Phys. Rev. B
93
,
165120
(
2016
).
9.
M.-C.
Kim
,
E.
Sim
, and
K.
Burke
,
Phys. Rev. Lett.
111
,
073003
(
2013
).
10.
X.
Zheng
,
M.
Liu
,
E. R.
Johnson
,
J.
Contreras-García
, and
W.
Yang
,
J. Chem. Phys.
137
,
214106
(
2012
).
11.
E. R.
Johnson
,
A.
Otero-de-la-Roza
, and
S. G.
Dale
,
J. Chem. Phys.
139
,
184116
(
2013
).
12.
D. J.
Tozer
and
F.
De Proft
,
J. Phys. Chem. A
109
,
8923
(
2005
).
13.
A. M.
Teale
,
F.
De Proft
, and
D. J.
Tozer
,
J. Chem. Phys.
129
,
044110
(
2008
).
14.
M. J. G.
Peach
,
A. M.
Teale
,
T.
Helgaker
, and
D. J.
Tozer
,
J. Chem. Theory Comput.
11
,
5262
(
2015
).
15.
A.
Ruzsinszky
,
J. P.
Perdew
,
G. I.
Csonka
,
O. A.
Vydrov
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
194112
(
2006
).
16.
A. D.
Dutoi
and
M.
Head-Gordon
,
Chem. Phys. Lett.
422
,
230
(
2006
).
17.
T.
Bally
and
G. N.
Sastry
,
J. Phys. Chem. A
101
,
7923
(
1997
).
18.
Y.
Zhang
and
W.
Yang
,
J. Chem. Phys.
109
,
2604
(
1998
).
19.
B. G.
Johnson
,
C. A.
Gonzales
,
P. M. W.
Gill
, and
J. A.
Pople
,
Chem. Phys. Lett.
221
,
100
(
1994
).
20.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
Phys. Rev. Lett.
100
,
146401
(
2008
).
21.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Phys. Rev. B
77
,
115123
(
2008
).
22.
H. J.
Kulik
,
M.
Cococcioni
,
D. A.
Scherlis
, and
N.
Marzari
,
Phys. Rev. Lett.
97
,
103001
(
2006
).
23.
G.
Ganzenmüller
,
N.
Berkaïne
,
A.
Fouqueau
,
M. E.
Casida
, and
M.
Reiher
,
J. Chem. Phys.
122
,
234321
(
2005
).
24.
A.
Droghetti
,
D.
Alfè
, and
S.
Sanvito
,
J. Chem. Phys.
137
,
124303
(
2012
).
25.
E. I.
Ioannidis
and
H. J.
Kulik
,
J. Chem. Phys.
143
,
034104
(
2015
).
26.
S. R.
Mortensen
and
K. P.
Kepp
,
J. Phys. Chem. A
119
,
4041
(
2015
).
27.
E. I.
Ioannidis
and
H. J.
Kulik
,
J. Phys. Chem. A
121
,
874
(
2017
).
28.
G.
Prokopiou
and
L.
Kronik
,
Chem. - Eur. J.
24
,
5173
(
2018
).
29.
P.
Verma
and
D. G.
Truhlar
,
J. Phys. Chem. Lett.
8
,
380
(
2017
).
30.
M. J.
Lucero
,
T. M.
Henderson
, and
G. E.
Scuseria
,
J. Phys.: Condens. Matter
24
,
145504
(
2012
).
31.
J.
Sun
,
A.
Ruzsinszky
, and
J. P.
Perdew
,
Phys. Rev. Lett.
115
,
036402
(
2015
).
32.
Y.
Zhao
and
D. G.
Truhlar
,
Theor. Chem. Acc.
120
,
215
(
2008
).
33.
N.
Mardirossian
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
16
,
9904
(
2014
).
34.
S.
Grimme
,
J. G.
Brandenburg
,
C.
Bannwarth
, and
A.
Hansen
,
J. Chem. Phys.
143
,
054107
(
2015
).
35.
B. G.
Janesko
,
E.
Proynov
,
J.
Kong
,
G.
Scalmani
, and
M. J.
Frisch
,
J. Phys. Chem. Lett.
8
,
4314
(
2017
).
36.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
37.
M. R.
Pederson
,
A.
Ruzsinszky
, and
J. P.
Perdew
, “
Communication: Self-interaction correction with unitary invariance in density functional theory
,”
J. Chem. Phys.
140
,
121103
(
2014
).
38.
M. R.
Pederson
,
J. Chem. Phys.
142
,
064112
(
2015
).
39.
K. P.
Withanage
,
K.
Trepte
,
J. E.
Peralta
,
T.
Baruah
,
R.
Zope
, and
K. A.
Jackson
,
J. Chem. Theory Comput.
14
,
4122
(
2018
).
40.
I.
Dabo
,
A.
Ferretti
,
N.
Poilvert
,
Y.
Li
,
N.
Marzari
, and
M.
Cococcioni
,
Phys. Rev. B
82
,
115121
(
2010
).
41.
S.
Lehtola
and
H.
Jónsson
,
J. Chem. Theory Comput.
10
,
5324
(
2014
).
42.
M.-C.
Kim
,
E.
Sim
, and
K.
Burke
,
J. Chem. Phys.
134
,
171103
(
2011
).
43.
M.-C.
Kim
,
H.
Park
,
S.
Son
,
E.
Sim
, and
K.
Burke
,
J. Phys. Chem. Lett.
6
,
3802
(
2015
).
44.
S.
Kümmel
and
L.
Kronik
,
Rev. Mod. Phys.
80
,
3
(
2008
).
45.
T.
Leininger
,
H.
Stoll
,
H.-J.
Werner
, and
A.
Savin
,
Chem. Phys. Lett.
275
,
151
(
1997
).
46.
J.
Toulouse
,
F.
Colonna
, and
A.
Savin
,
Phys. Rev. A
70
,
062505
(
2004
).
47.
O. A.
Vydrov
and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
234109
(
2006
).
48.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
J. Chem. Phys.
126
,
191109
(
2007
).
49.
H.
Iikura
,
T.
Tsuneda
,
T.
Yanai
, and
K.
Hirao
,
J. Chem. Phys.
115
,
3540
(
2001
).
50.
R.
Baer
and
D.
Neuhauser
,
Phys. Rev. Lett.
94
,
043002
(
2005
).
51.
T.
Tsuneda
,
J.-W.
Song
,
S.
Suzuki
, and
K.
Hirao
,
J. Chem. Phys.
133
,
174101
(
2010
).
52.
T. M.
Henderson
,
B. G.
Janesko
, and
G. E.
Scuseria
,
J. Phys. Chem. A
112
,
12530
(
2008
).
53.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
54.
A. V.
Krukau
,
O. A.
Vydrov
,
A. F.
Izmaylov
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
224106
(
2006
).
55.
M. A.
Rohrdanz
,
K. M.
Martins
, and
J. M.
Herbert
,
J. Chem. Phys.
130
,
054112
(
2009
).
56.
V. I.
Anisimov
,
J.
Zaanen
, and
O. K.
Andersen
,
Phys. Rev. B
44
,
943
(
1991
).
57.
A. I.
Liechtenstein
,
V. I.
Anisimov
, and
J.
Zaanen
,
Phys. Rev. B
52
,
R5467
(
1995
).
58.
S. L.
Dudarev
,
G. A.
Botton
,
S. Y.
Savrasov
,
C. J.
Humphreys
, and
A. P.
Sutton
,
Phys. Rev. B
57
,
1505
(
1998
).
59.
H. J.
Kulik
,
J. Chem. Phys.
142
,
240901
(
2015
).
60.
M.
Cococcioni
and
S.
de Gironcoli
,
Phys. Rev. B
71
,
035105
(
2005
).
61.
T. Z. H.
Gani
and
H. J.
Kulik
,
J. Chem. Theory Comput.
13
,
5443
(
2017
).
62.
E.
Livshits
and
R.
Baer
,
Phys. Chem. Chem. Phys.
9
,
2932
(
2007
).
63.
L.
Kronik
,
T.
Stein
,
S.
Refaely-Abramson
, and
R.
Baer
,
J. Chem. Theory Comput.
8
,
1515
(
2012
).
64.
M.
Srebro
and
J.
Autschbach
,
J. Phys. Chem. Lett.
3
,
576
(
2012
).
65.
T.
Körzdörfer
and
J.-L.
Brédas
,
Acc. Chem. Res.
47
,
3284
(
2014
).
66.
J.
Autschbach
and
M.
Srebro
,
Acc. Chem. Res.
47
,
2592
(
2014
).
67.
J. D.
Gledhill
,
M. J. G.
Peach
, and
D. J.
Tozer
,
J. Chem. Theory Comput.
9
,
4414
(
2013
).
68.
U.
Salzner
and
R.
Baer
,
J. Chem. Phys.
131
,
231101
(
2009
).
69.
Q.
Zhao
,
E. I.
Ioannidis
, and
H. J.
Kulik
,
J. Chem. Phys.
145
,
054109
(
2016
).
70.
A.
Mahler
,
B. G.
Janesko
,
S.
Moncho
, and
E. N.
Brothers
,
J. Chem. Phys.
148
,
244106
(
2018
).
71.
A.
Bajaj
,
J. P.
Janet
, and
H. J.
Kulik
,
J. Chem. Phys.
147
,
191101
(
2017
).
72.
N. Q.
Su
,
C.
Li
, and
W.
Yang
,
Proc. Natl. Acad. Sci. U. S. A.
115
,
9678
(
2018
).
73.
L. A.
Agapito
,
S.
Curtarolo
, and
M. B.
Nardelli
,
Phys. Rev. X
5
,
011006
(
2015
).
74.
T. Z.
Gani
and
H. J.
Kulik
,
J. Chem. Theory Comput.
12
,
5931
(
2016
).
75.
Q.
Zhao
and
H.
Kulik
,
J. Chem. Theory Comput.
14
,
670
(
2018
).
76.
T. J.
Duignan
and
J.
Autschbach
,
J. Chem. Theory Comput.
12
,
3109
(
2016
).
77.
A. D.
Becke
,
J. Chem. Phys.
138
,
074109
(
2013
).
78.
J. P.
Perdew
,
R. G.
Parr
,
M.
Levy
, and
J. L.
Balduz
,
Phys. Rev. Lett.
49
,
1691
(
1982
).
79.
J. P.
Perdew
and
M.
Levy
,
Phys. Rev. Lett.
51
,
1884
(
1983
).
80.
L. J.
Sham
and
M.
Schlüter
,
Phys. Rev. Lett.
51
,
1888
(
1983
).
81.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
Phys. Rev. Lett.
102
,
066403
(
2009
).
82.
A. J.
Cohen
,
P.
Mori-Sanchez
, and
W. T.
Yang
,
Chem. Rev.
112
,
289
(
2012
).
83.
J. J.
Phillips
,
A. A.
Kananenka
, and
D.
Zgid
,
J. Chem. Phys.
142
,
194108
(
2015
).
84.
J. T.
Margraf
and
R.
Bartlett
,
J. Chem. Phys.
148
,
221103
(
2018
).
85.
E. R.
Johnson
and
J.
Contreras-García
,
J. Chem. Phys.
135
,
081103
(
2011
).
86.
E.
Sagvolden
and
J. P.
Perdew
,
Phys. Rev. A
77
,
012517
(
2008
).
87.
P.
Mori-Sanchez
and
A. J.
Cohen
,
Phys. Chem. Chem. Phys.
16
,
14378
(
2014
).
88.
A.
Kramida
,
Y.
Ralchenko
,
J.
Reader
, and
N. A.
Team
, NIST Atomic Spectra Database (verision 5.5.6), [Online]. Available: https://physics.nist.gov/asd.
89.
X.
Yang
,
A.
Patel
,
R.
Miranda-Quintana
,
F.
Heidar-Zadeh
,
C.
Gonzalez-Espinoza
, and
P.
Ayers
,
J. Chem. Phys.
145
,
031102
(
2016
).
90.
U.
von Barth
,
Phys. Rev. A
20
,
1693
(
1979
).
91.
A.
Cohen
,
P.
Mori-Sanchez
, and
W.
Yang
,
J. Chem. Phys.
129
,
121104
(
2008
).
92.
B.
Himmetoglu
,
A.
Floris
,
S.
De Gironcoli
, and
M.
Cococcioni
,
Int. J. Quantum Chem.
114
,
14
(
2014
).
93.
T.
Stein
,
J.
Autschbach
,
N.
Govind
,
L.
Kronik
, and
R.
Baer
,
J. Phys. Chem. Lett.
3
,
3740
(
2012
).
94.
H. J.
Kulik
and
N.
Marzari
,
J. Chem. Phys.
129
,
134314
(
2008
).
95.
W.
Pickett
,
S.
Erwin
, and
E.
Ethridge
,
Phys. Rev. B
58
,
1201
(
1998
).
96.
P.
Giannozzi
,
O.
Andreussi
,
T.
Brumme
,
O.
Bunau
,
M. B.
Nardelli
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
M.
Cococcioni
,
N.
Colonna
,
I.
Carnimeo
,
A.
Dal Corso
,
S.
de Gironcoli
,
P.
Delugas
,
R. A.
DiStasio
, Jr.
,
A.
Ferretti
,
A.
Floris
,
G.
Fratesi
,
G.
Fugallo
,
R.
Gebauer
,
U.
Gerstmann
,
F.
Giustino
,
T.
Gorni
,
J.
Jia
,
M.
Kawamura
,
H.-Y.
Ko
,
A.
Kokalj
,
E.
Küçükbenli
,
M.
Lazzeri
,
M.
Marsili
,
N.
Marzari
,
F.
Mauri
,
N. L.
Nguyen
,
H.-V.
Nguyen
,
A.
Otero-de-la-Roza
,
L.
Paulatto
,
S.
Poncé
,
D.
Rocca
,
R.
Sabatini
,
B.
Santra
,
M.
Schlipf
,
A. P.
Seitsonen
,
A.
Smogunov
,
I.
Timrov
,
T.
Thonhauser
,
P.
Umari
,
N.
Vast
,
X.
Wu
, and
S.
Baroni
,
J. Phys.: Condens. Matter
29
,
465901
(
2017
).
97.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
98.
D.
Vanderbilt
,
Phys. Rev. B
41
,
7892
(
1990
).
99.
D. R.
Hamann
,
M.
Schlüter
, and
C.
Chiang
,
Phys. Rev. Lett.
43
,
1494
(
1979
).
100.

Quantum-ESPRESSO, Quantum-ESPRESSO Pseudopotentials.

101.
G. J.
Martyna
and
M. E.
Tuckerman
,
J. Chem. Phys.
110
,
2810
(
1999
).
102.
MATLAB version 9.4.0 (
2018
).
103.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
104.
A. I.
Krylov
,
C. D.
Sherrill
, and
M.
Head-Gordon
,
J. Chem. Phys.
113
,
6509
(
2000
).
105.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
106.
W.
Yang
,
Y.
Zhang
, and
P. W.
Ayers
,
Phys. Rev. Lett.
84
,
5172
(
2000
).
107.
D.
Hait
and
M.
Head-Gordon
,
J. Phys. Chem. Lett.
9
,
6280
(
2018
).
109.
M. G.
Medvedev
,
I. S.
Bushmarinov
,
J.
Sun
,
J. P.
Perdew
, and
K. A.
Lyssenko
,
Science
355
,
49
(
2017
).
110.
K. R.
Brorsen
,
Y.
Yang
,
M. V.
Pak
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. Lett.
8
,
2076
(
2017
).
111.
P. D.
Mezei
,
G. I.
Csonka
, and
M.
Kállay
,
J. Chem. Theory Comput.
13
,
4753
(
2017
).
112.
P.
Verma
,
A.
Perera
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
524
,
10
(
2012
).
113.
I. E.
Brumboiu
,
G.
Prokopiou
,
L.
Kronik
, and
B.
Brena
,
J. Chem. Phys.
147
,
044301
(
2017
).
114.
T. J.
Duignan
,
J.
Autschbach
,
E.
Batista
, and
P.
Yang
,
J. Chem. Theory Comput.
13
,
3614
(
2017
).
115.
F.
Jensen
,
J. Chem. Theory Comput.
6
,
2726
(
2010
).
116.
D. N.
Komsa
and
V. N.
Staroverov
,
J. Chem. Theory Comput.
12
,
5361
(
2016
).

Supplementary Material