The flat-plane condition is the union of two exact constraints in electronic structure theory: (i) energetic piecewise linearity with fractional electron removal or addition and (ii) invariant energetics with change in electron spin in a half filled orbital. Semi-local density functional theory (DFT) fails to recover the flat plane, exhibiting convex fractional charge errors (FCE) and concave fractional spin errors (FSE) that are related to delocalization and static correlation errors. We previously showed that DFT+U eliminates FCE but now demonstrate that, like other widely employed corrections (i.e., Hartree-Fock exchange), it worsens FSE. To find an alternative strategy, we examine the shape of semi-local DFT deviations from the exact flat plane and we find this shape to be remarkably consistent across ions and molecules. We introduce the judiciously modified DFT (jmDFT) approach, wherein corrections are constructed from few-parameter, low-order functional forms that fit the shape of semi-local DFT errors. We select one such physically intuitive form and incorporate it self-consistently to correct semi-local DFT. We demonstrate on model systems that jmDFT represents the first easy-to-implement, no-overhead approach to recovering the flat plane from semi-local DFT.

Approximate density functional theory (DFT) is one of the most widely used electronic structure methods, despite the fact that presently available pure exchange-correlation (xc) approximations in Kohn-Sham DFT are plagued by both one- and many-electron self-interaction errors (SIE),1–5 which give rise to well-known problems in dissociation energies,2,6–9 barrier heights,10 band gaps,11,12 electron affinities,13–15 and other manifestations of delocalization error.16–18 A number of generalizations to Kohn-Sham DFT, e.g., range-separated hybrids,19–25 global hybrids, and DFT+U,26–29 mitigate errors in part by recovering the piecewise linearity of the energy with respect to fractional electron removal or addition, a requirement of any exact electronic structure method.30–32 The combination of piecewise linearity with the additional constraint that the energy is invariant with respect to the spin of an electron in isoenergetic orbitals produces the flat-plane condition33,34 in which two planes meet at a fractional spin line35 (FSL) seam.

The flat-plane condition applies to any three electronic states with N − 1, N, and N + 1 electrons, where (i) an orbital is half-occupied, doubly occupied, and empty in the N-,N + 1-, and N − 1-electron states and (ii) the N-electron multiplicity is one higher (e.g., doublet) than the N − 1/N + 1-electron multiplicity (e.g., singlet). Ionization potentials (IPs) or electron affinities (EAs) connecting these electronic states, which are available for ions from experimental databases,36 enable construction of “exact” flat planes. An exhaustive study37 of such flat planes has revealed that nearly all feasible electronic state combinations preserve this shape, with a few exceptions for highly ionized and excited electronic states outside the scope of the current work. Our focus on ground state DFT narrows the study to cases where all three electronic states correspond to the ground state of the ion or to the lowest state with that quantum number: ionization of (i) s electrons in H–Rb or He+–Sr+ (N − 1/N + 1: 1S, N: 2S), (ii) p electrons in N–Bi or O+–Te+ (N − 1/N + 1: 3P, N: 4S), (iii) d electrons in Fe3+ or Co4+ (N − 1/N + 1: 5D, N: 6S) (Text S1 and Table S1 of the supplementary material). We take a similar approach to construct molecular (H2, Na2) flat planes from correlated wavefunction theory reference IP and EA data (Texts S1 and S2 of the supplementary material).

We recently demonstrated38 that DFT+U26–29 recovers piecewise linearity30–32 and the derivative discontinuity30 in representative transition metal ions and inorganic complexes. The most widely applied DFT+U correction26–29,39 is

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

where UnlI is applied to the nl subshell of the Ith atom, σ is the spin index, and nnlI,σ is the occupation matrix with elements, nmmI,σ=vψvϕmIϕmIψv, obtained from projecting all extended states (i.e., with index ν), ψ, onto atomic orbitals, ϕ. DFT+U thus acts as a local correction on valence orbitals with strong nl character. In Eq. (1), differences in Coulomb repulsion between the same- and opposite-spin electrons (i.e., exchange) have been neglected, and UnlI is an effective Ueff = UJ.28,29,39 Even in the atomic limit where electron count is reflected in occupations (i.e., δN ∼ δn), piecewise linearity is not recovered at linear-response Ueff,39,40 i.e., 2En2, but is recovered when the average curvature from frontier orbitals,41 i.e., 2EN2, is applied.

By construction, this DFT+U formulation should dramatically increase energetic errors along the FSL (i.e., nα + nβ = 1, where nαand nβ are the fractional electron counts of α and β orbitals emptying or filling in the flat plane). As an example, recovering piecewise linearity with DFT+U with fractional electron addition from He+ to He [i.e., the fractional charge line (FCL)] requires Ueff = 2EN2 = 16 eV (Fig. 1). Examining errors along the He+ FSL reveals that modest errors for semi-local DFT (here, PBE52 max. value ca. 2 eV) are dramatically (i.e., 200%) increased to 6 eV with DFT+U (Fig. 1). Similar increases in FSL errors are known to occur for methods that incorporate Hartree-Fock (HF) exchange, e.g., full HF, global hybrids, or range-separated hybrids,35,42 motivating alternative strategies.43 Such a large increase in fractional spin error has been largely overlooked in the application of tuned range-separated hybrids for organic molecules44–49 and some inorganic50,51 systems, with few exceptions.42 It has been known for some time that increasing FSL error is directly related to increasing static correlation errors (SCEs).35 This analysis highlights that, like incorporation of HF exchange, DFT+U directly increases SCE, which should be particularly concerning for already-SCE-prone transition metal complexes and solids to which DFT+U is widely applied.

FIG. 1.

Deviation from linearity (left) and relative FSL energy (right) of PBE (gray squares) and with Ueff = 16 eV (blue circles) for the He+ to He line. The maximum deviation of PBE is annotated. Lines correspond to interpolation of the explicitly calculated points.

FIG. 1.

Deviation from linearity (left) and relative FSL energy (right) of PBE (gray squares) and with Ueff = 16 eV (blue circles) for the He+ to He line. The maximum deviation of PBE is annotated. Lines correspond to interpolation of the explicitly calculated points.

Close modal

Thus, we depart from traditional methods to identify if other functional forms can recover piecewise linearity and the derivative discontinuity while decreasing or even eliminating FSL errors. We consider functional forms that are easy to implement as corrections to the total energy and the potential within the context of a Kohn-Sham DFT calculation without increasing computational cost. In order to identify such functional forms, we return to the flat-plane condition and evaluate the deviations of semi-local DFT functionals from the exact flat-plane condition for ions and molecules (here, PBE52; behavior from other functionals shown in Fig. S1 of the supplementary material). We choose He+ as our primary worked example as observations are transferable to other elements but the large PBE errors make it a challenging test case (Fig. 2 and Texts S1–S4 and Tables S1–S14 of the supplementary material).

FIG. 2.

Flat-plane error of the PBE functional (in eV) compared to experiment (schematics of PBE and experiment shown in inset) colored according to sign of errors (from negative in red through white to positive in blue). Energies are aligned at He+.

FIG. 2.

Flat-plane error of the PBE functional (in eV) compared to experiment (schematics of PBE and experiment shown in inset) colored according to sign of errors (from negative in red through white to positive in blue). Energies are aligned at He+.

Close modal

We align the N-electron (He+) PBE and exact energies, distributing IP/EA errors into the N + 1 and N − 1 endpoints (errors in Table S1 of the supplementary material). The errors in PBE with respect to the exact flat plane are convex along the FCL and concave along the FSL (Fig. 2 and Fig. S2 of the supplementary material). Examination of this error also (i) suggests that the smoothly varying deviations are likely to be suitably fit with a low-order polynomial of nα and nβ in the spirit of DFT+U and (ii) highlights why a penalty on fractional occupations [Eq. (1)] can correct the edges of the flat plane (i.e., the FCLs) but will dome the middle region, exacerbating FSL errors (Fig. 2). Advantages of low-order polynomial corrections are that (i) they add no overhead in the solid state and (ii) they are transparent, enabling the full derivation of their behavior over a range of expected electron configurations,38 as opposed to more opaque many-parameter approximate functionals that can unpredictably lack good extrapolative performance on molecules that differ from training data.

First focusing on the FSL error, it is clear that a second-order expression is needed that (i) opposes the intrinsic PBE penalty for nα = 1/2, nβ = 1/2 with respect to an nα = 1 or nβ = 1 configuration and (ii) approaches zero as nα or nβ approach 1 (Figs. 1 and 2). A functional form that can eliminate this FSL error is the J correction in DFT+U+J,

EDFT+U+J=EDFT+12I,σnl(UnlIJnlI)[Tr(nnlI,σ(1nnlI,σ))]+12I,σnlJnlI[Tr(nnlI,σnnlI,σ)].
(2)

If we apply Eq. (2) to only correct the FSL in He+, we would thus require a negative J of −16 eV and Ueff = UJ = 0. There is a key reason why DFT+U+J alone cannot recover the flat-plane condition that should be relevant to researchers wishing to broadly employ the method for correlated molecules: the J term in Eq. (2) increases with increasing filling of the nl subshell, producing an arbitrary shift in the energy of differing oxidation states (Fig. S3 of the supplementary material). The primary assumption in most functional tuning approaches, consistent with our own data, is that integer-electron endpoint errors are modest and endpoint adjusting terms (e.g., J) are not beneficial in most cases where the magnitude of the needed correction is not known a priori (Table S1 of the supplementary material).

Thus, we now consider generalizations of DFT+U+J that can be used to construct corrections to semi-local functionals on the basis of the shape of the error, much as a jello mold is used to create a 3D shape. The resulting form can then be incorporated into a self-consistent field calculation to counteract such errors and recover the flat plane. Working with ions and molecules in equilibrium (e.g., H2, Na2), we use as the input to our fit the difference between semi-local DFT and exact flat-plane energies along a grid of nα and nβ values for the emptying or filling electron level. We judge these corrections by how well they fit to the shape of the error, through root mean square error (RMSE) measures and by the number of new or established (e.g., U, J) parameters they introduce. We assume changes in electron count will be reflected in occupations. Low order polynomials will require two sets of coefficients on either side of the plane to introduce a discontinuity. Our approach, which we refer to in its broadest form as judiciously modified DFT (jmDFT), no longer requires a grounding in the Hubbard model in its original form53 or extensions54–58 to construct model potentials. Consistent with earlier observations, DFT+U or DFT+U+J terms make poor fits to the PBE errors: the best fits are to DFT+U+J with differing U and J values (i.e., 4 parameters) at around 0.1-0.5 eV RMSE, and DFT+U-only errors are much larger (Tables S3–S14 of the supplementary material).

Examining shapes and parameters of fit expressions across elements confirms that errors are generally the same shape, but this also leads to observations consistent with earlier work on only delocalization error:38,45,59–62 heavier elements exhibit reduced convexity and concavity as the addition or removal of a single electron becomes a smaller fraction of total electron count (i.e., smaller coefficients in fits, Tables S3–S14 of the supplementary material).

The most general fit can obtain an order of magnitude lower RMSE (ca. 0.01-0.05 eV overall, 0.08 eV in He, Fig. 3) than DFT+U-based approaches by applying a second-order polynomial in the electron count of orbitals that empty and fill between states in the flat plane, nα and nβ,

ΔE=a+bv1122+cv2+dv22=a+b4(nαnβ)2+c2nα+nβ1+d4nα+nβ12.
(3)

This functional form arises from observing that transforming the flat-plane error to a coordinate system, v1 = 1/2(nα − nβ + 1) and v2 = 1/2(nα + nβ − 1), aligns the FSL to v2 = 0 and reveals the directions of approximately parabolic dependence (Text S3 of the supplementary material). Rather than fitting the three extrema with a fourth order polynomial with up to fifteen coefficients, we employ asymmetric coefficients (i.e., 8 parameters) in Eq. (3) to introduce a discontinuity at the FSL (Text S3 of the supplementary material). Although this leads to a good fit, the large number and non-intuitive nature of these parameters, including linear terms63 and constants, leads us to disfavor this correction form. This most general form bears no connection to DFT+U+J, but grouping terms reveals U- or J-like components (Text S3 of the supplementary material).

FIG. 3.

(a) Needed correction (ΔPBE, in eV) to PBE with dots colored as in the legend in (b) along with the U+J/J′ correction in the translucent blue surface. (b) Cumulative absolute error (eV) over the plane for U+J, U+J/J′, U+J+K, and polynomial (poly). The stacked bar is annotated by fit RMSE, and components are (i) He+ to He FCL (FCL0, red), (ii) the nα + nβ > 1 plane (orange), (iii) He2+ to He+ FCL (FCL+, blue), (iv) the nα + nβ < 1 plane (green), and (v) FSL (gray).

FIG. 3.

(a) Needed correction (ΔPBE, in eV) to PBE with dots colored as in the legend in (b) along with the U+J/J′ correction in the translucent blue surface. (b) Cumulative absolute error (eV) over the plane for U+J, U+J/J′, U+J+K, and polynomial (poly). The stacked bar is annotated by fit RMSE, and components are (i) He+ to He FCL (FCL0, red), (ii) the nα + nβ > 1 plane (orange), (iii) He2+ to He+ FCL (FCL+, blue), (iv) the nα + nβ < 1 plane (green), and (v) FSL (gray).

Close modal

We return to DFT+U+J to identify additions and adjustments that can recover the flat plane. We first add a third term, K, that counteracts energy increases from J with orbital filling. The K fitting coefficients will be zero for the N − 1- to N-electron portion of the flat plane but counteract the J term on the N- to N + 1-electron side (He2+ to He+: K = 0.6 eV, He+ to He: K = −25.5 eV),

EK=12I,σnlKnlITr(nnlI,σnnlI,σ+1)(nnlI,σnnlI,σ+1).
(4)

Physically, Eq. (4) could provide a finer tuning of spin state ordering than DFT+U alone by having no effect on intermediate spin states but favoring low- or high-spin states, if K is negative or positive, respectively. Asymmetric DFT+U+J+K introduces six parameters and relies on the difference of two large effects (i.e., +J and +K) to cancel, introducing precision or stability issues but otherwise giving a more intuitive fit with comparable quality to the polynomial (Fig. 3, asym. RMSE in He: 0.11 eV, 0.01-0.05 eV for heavier elements, see Tables S3–S14 and Fig. S4 of the supplementary material). Alternatively, we adjust the J term in Eq. (2) (U+J/J′), resetting the electron count on the N- to N + 1-electron side,

EJ=12I,σnlJnlITrnnlI,σ1nnlI,σ1,
(5)

retaining Eq. (2) for the N − 1- to N-electron side. Although the difference in curvature on the FCL is known to typically require different parameters for different ionization processes,38,64 the U+J/J′ form appears to permit reduction in the number of parameters: symmetric coefficients (i.e., J = J′, 2 parameters) are of comparable quality to asymmetric (4 parameters, see Tables S3–S14 of the supplementary material). The symmetric coefficient model provides among the lowest RMSEs (ca. 0.03 eV in most s-block elements and 0.10 eV in most p-block elements, higher in He: sym. RMSE = 0.35 eV, asym. 0.15 eV) (Fig. 3 and see Tables S3–S14 of the supplementary material). Here, U is positive, whereas J is negative, and increasing atomic number decreases coefficients (i.e., for He sym. U+J/J′ U = 20.4 and J = −30.4 eV and for Ca sym. U+J/J′ U = 4.8 and J = −5.7 eV, see Table S5 of the supplementary material). Alternative functional forms considered were not pursued further because they introduced more parameters or did not improve upon standard DFT+U+J (Text S3 and Tables S3–S14 of the supplementary material).

Decomposing total errors into contributions from different regions quantitatively confirms earlier observations on four representative fits for He (Fig. 3). The lowest errors observed in He+ with the polynomial are achieved through balanced error elimination in all regions of the plane but especially eliminating FSL errors (RMSE = 0.08 eV). Conversely, U+J produces large errors there and in the He+ to He FCL (RMSE = 0.56 eV) that are eliminated with U+J/J′ (RMSE = 0.15 eV, Fig. 3). For the He example, the U+J/J′ RMSE (0.15 eV) remains larger than the 6-parameter U+J+K (RMSE = 0.11 eV) and polynomial (RMSE = 0.08 eV) correction because these approaches correct integer electron errors, improving the He2+ to He+ FCL [PBE He+ IP error of 0.85 eV, Fig. 3(a)]. The three forms are more comparable in other elements that have smaller endpoint errors (Tables S3–S14 of the supplementary material).

At this point, we select the U+J/J′ correction as the functional form that we will employ in our jmDFT approach to build a fully self-consistent functional that can be applied across the periodic table. We chose this form because in several cases (e.g., s-block Be, Mg, Ca or p-block S, Se, P, As), the symmetric, 2-parameter fit is of comparable fit to the asymmetric, 4-parameter one and it does not introduce an adjustable parameter to the total energy at integer electron number (Tables S3–S14 of the supplementary material). The energy expression in Eqs. (2) and (5) is incorporated self-consistently by inclusion of a potential correction of the form

VjmDFT=I,σnlmUnlI212nnl,mI,σ+Jnnl,mI,σ1ϕnl,mIϕnl,mI,
(6)

where 1 in the second parentheses is only present for the J′ case. We have added both this correction and the U+J+K form into locally modified copies of Quantum-ESPRESSO,65 the source code available upon request. Incorporating this new potential and energy correction introduces no overhead or complexity beyond the ingredients of widely employed semi-local DFT corrections, i.e., we require that occupations obtained from projecting onto atomic states remain a good proxy for the fractional electron count that we used in the initial fit (i.e., δN ∼ δn). We recently demonstrated the range of applicability of this assumption in larger complexes, and alternate forms, currently being developed by our group, could be employed more suitably in the rare cases where this approximation fails.38 Self-consistent evaluation of the energy with the jmDFT potential correction (N − 1 to N: U = −12, J = −36, i.e., Ueff = 24 eV; N to N + 1: U = −8, J′ = −25, i.e., Ueff = 17 eV) satisfied this requirement at all fractional charge and spin points for our He example (other examples in Figs. S5–S9 of the supplementary material). The final plot of the energy versus nα and nβ electron filling reveals a recovered flat plane (Fig. 4). This flat plane includes the seam at the FSL and a derivative discontinuity at the meeting point between the two FCLs (Fig. 4). Examining the FSL and FCL from He+ to He reveals that errors obtained with PBE are dramatically reduced with respect to the exact answer for He. Evaluation of a number of other s electron systems reveals the same result, with the only difference being the magnitude of the flat-plane deviations to correct, which generally decreases with increasing electron number or in stretched molecular geometries (e.g., stretched H2, Figs. S5–S9 and Table S14 of the supplementary material). Higher angular momenta (e.g., 2p) introduce additional conditions,35,66 but jmDFT should be transferable (Fig. S10 of the supplementary material).

FIG. 4.

(Left) jmDFT energies annotated (1, FSL and 2, FCL0) for insets on the right. (Right) The FSL (top) and FCL0 (bottom) for jmDFT (green circles), PBE (gray squares), and exact result (open circles), with the legend shown on the left.

FIG. 4.

(Left) jmDFT energies annotated (1, FSL and 2, FCL0) for insets on the right. (Right) The FSL (top) and FCL0 (bottom) for jmDFT (green circles), PBE (gray squares), and exact result (open circles), with the legend shown on the left.

Close modal

In conclusion, we introduce a practical, no-overhead solution to recovering the flat plane within semi-local DFT. This jmDFT framework involves fitting to known errors with electron count and incorporating the correction into the Kohn-Sham potential at no additional computational cost. We have identified that appropriate functional forms are remarkably transferable across ions and small molecules, as judged through similarly low RMSEs for expressions that capture FCL convexity and FSL concavity. This result suggests that one general functional form should be suitable to recover exact conditions in a wide range of small molecules. We selected a form reminiscent of Hubbard model corrections because it balances the fewest parameters with the physical transparency in the correction. As we extend our methodology to larger systems, stretched geometries, or higher angular momenta, we anticipate that components in this methodology, such as the projections and the polynomial correction, may be different. We expect that these functional forms may help guide alternate approaches to SCE minimization in DFT.

See supplementary material for flat-plane construction and simulation details, experimental reference data, comparison of other xcs, list of pseudopotentials, other fitting expressions and procedure, RMSE and fitting coefficients for six models and s-, p-, and d-block elements, example PBE errors for representative systems, components of correction terms, recovery of the flat plane in three other atoms and an equilibrium geometry molecule, and extensions to higher angular momenta or stretched molecules.

The authors acknowledge support by the Department of Energy under Grant No. DE-SC0018096, the Office of Naval Research under Grant No. N00014-17-1-2956, and MIT Energy Initiative seed grants (2014, 2017). J.P.J. was supported in part by an MIT-SUTD Graduate fellowship. H.J.K. holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund. This work was carried out in part using computational resources from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575. The authors thank Adam H. Steeves for providing a critical reading of the manuscript.

1.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
J. Chem. Phys.
125
,
201102
(
2006
).
2.
A.
Ruzsinszky
,
J. P.
Perdew
,
G. I.
Csonka
,
O. A.
Vydrov
, and
G. E.
Scuseria
,
J. Chem. Phys.
126
,
104102
(
2007
).
3.
R.
Haunschild
,
T. M.
Henderson
,
C. A.
Jiménez-Hoyos
, and
G. E.
Scuseria
,
J. Chem. Phys.
133
,
134116
(
2010
).
4.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Science
321
,
792
(
2008
).
5.
T.
Schmidt
and
S.
Kümmel
,
Phys. Rev. B
93
,
165120
(
2016
).
6.
A.
Ruzsinszky
,
J. P.
Perdew
,
G. I.
Csonka
,
O. A.
Vydrov
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
194112
(
2006
).
7.
A. D.
Dutoi
and
M.
Head-Gordon
,
Chem. Phys. Lett.
422
,
230
(
2006
).
8.
T.
Bally
and
G. N.
Sastry
,
J. Phys. Chem. A
101
,
7923
(
1997
).
9.
Y.
Zhang
and
W.
Yang
,
J. Chem. Phys.
109
,
2604
(
1998
).
10.
B. G.
Johnson
,
C. A.
Gonzales
,
P. M. W.
Gill
, and
J. A.
Pople
,
Chem. Phys. Lett.
221
,
100
(
1994
).
11.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
Phys. Rev. Lett.
100
,
146401
(
2008
).
12.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Phys. Rev. B
77
,
115123
(
2008
).
13.
D. J.
Tozer
and
F.
De Proft
,
J. Phys. Chem. A
109
,
8923
(
2005
).
14.
A. M.
Teale
,
F.
De Proft
, and
D. J.
Tozer
,
J. Chem. Phys.
129
,
044110
(
2008
).
15.
M. J. G.
Peach
,
A. M.
Teale
,
T.
Helgaker
, and
D. J.
Tozer
,
J. Chem. Theory Comput.
11
,
5262
(
2015
).
16.
M.-C.
Kim
,
E.
Sim
, and
K.
Burke
,
Phys. Rev. Lett.
111
,
073003
(
2013
).
17.
X.
Zheng
,
M.
Liu
,
E. R.
Johnson
,
J.
Contreras-García
, and
W.
Yang
,
J. Chem. Phys.
137
,
214106
(
2012
).
18.
E. R.
Johnson
,
A.
Otero-de-la-Roza
, and
S. G.
Dale
,
J. Chem. Phys.
139
,
184116
(
2013
).
19.
T.
Leininger
,
H.
Stoll
,
H.-J.
Werner
, and
A.
Savin
,
Chem. Phys. Lett.
275
,
151
(
1997
).
20.
J.
Toulouse
,
F.
Colonna
, and
A.
Savin
,
Phys. Rev. A
70
,
062505
(
2004
).
21.
O. A.
Vydrov
and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
234109
(
2006
).
22.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
J. Chem. Phys.
126
,
191109
(
2007
).
23.
H.
Iikura
,
T.
Tsuneda
,
T.
Yanai
, and
K.
Hirao
,
J. Chem. Phys.
115
,
3540
(
2001
).
24.
R.
Baer
and
D.
Neuhauser
,
Phys. Rev. Lett.
94
,
043002
(
2005
).
25.
T.
Tsuneda
,
J.-W.
Song
,
S.
Suzuki
, and
K.
Hirao
,
J. Chem. Phys.
133
,
174101
(
2010
).
26.
V. I.
Anisimov
,
J.
Zaanen
, and
O. K.
Andersen
,
Phys. Rev. B
44
,
943
(
1991
).
27.
A. I.
Liechtenstein
,
V. I.
Anisimov
, and
J.
Zaanen
,
Phys. Rev. B
52
,
R5467
(
1995
).
28.
S. L.
Dudarev
,
G. A.
Botton
,
S. Y.
Savrasov
,
C. J.
Humphreys
, and
A. P.
Sutton
,
Phys. Rev. B
57
,
1505
(
1998
).
29.
H. J.
Kulik
,
J. Chem. Phys.
142
,
240901
(
2015
).
30.
J. P.
Perdew
,
R. G.
Parr
,
M.
Levy
, and
J. L.
Balduz
,
Phys. Rev. Lett.
49
,
1691
(
1982
).
31.
J. P.
Perdew
and
M.
Levy
,
Phys. Rev. Lett.
51
,
1884
(
1983
).
32.
L. J.
Sham
and
M.
Schlüter
,
Phys. Rev. Lett.
51
,
1888
(
1983
).
33.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
Phys. Rev. Lett.
102
,
066403
(
2009
).
34.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Chem. Rev.
112
,
289
(
2011
).
35.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
J. Chem. Phys.
129
,
121104
(
2008
).
36.
A.
Kramida
,
Ralchenko
,
Yu.
,
Reader
,
J.
and
NIST ASD Team
, NIST Atomic Spectra Database (version 5.3) (Online),
2015
.
37.
X. D.
Yang
,
A. H.
Patel
,
R. A.
Miranda-Quintana
,
F.
Heidar-Zadeh
,
C. E.
González-Espinoza
, and
P. W.
Ayers
,
J. Chem. Phys.
145
,
031102
(
2016
).
38.
Q.
Zhao
,
E. I.
Ioannidis
, and
H. J.
Kulik
,
J. Chem. Phys.
145
,
054109
(
2016
).
39.
M.
Cococcioni
and
S.
de Gironcoli
,
Phys. Rev. B
71
,
035105
(
2005
).
40.
H. J.
Kulik
,
M.
Cococcioni
,
D. A.
Scherlis
, and
N.
Marzari
,
Phys. Rev. Lett.
97
,
103001
(
2006
).
41.
T.
Stein
,
J.
Autschbach
,
N.
Govind
,
L.
Kronik
, and
R.
Baer
,
J. Phys. Chem. Lett.
3
,
3740
(
2012
).
42.
T. J.
Duignan
and
J.
Autschbach
,
J. Chem. Theory Comput.
12
,
3109
(
2016
).
43.
E. R.
Johnson
and
J.
Contreras-García
,
J. Chem. Phys.
135
,
081103
(
2011
).
44.
I.
Tamblyn
,
S.
Refaely-Abramson
,
J. B.
Neaton
, and
L.
Kronik
,
J. Phys. Chem. Lett.
5
,
2734
(
2014
).
45.
A.
Karolewski
,
L.
Kronik
, and
S.
Kümmel
,
J. Chem. Phys.
138
,
204115
(
2013
).
46.
J.
Autschbach
and
M.
Srebro
,
Acc. Chem. Res.
47
,
2592
(
2014
).
47.
T.
Körzdörfer
and
J.-L.
Brédas
,
Acc. Chem. Res.
47
,
3284
(
2014
).
48.
S.
Refaely-Abramson
,
R.
Baer
, and
L.
Kronik
,
Phys. Rev. B
84
,
075144
(
2011
).
49.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
N.
Govind
,
J.
Autschbach
,
J. B.
Neaton
,
R.
Baer
, and
L.
Kronik
,
Phys. Rev. Lett.
109
,
226405
(
2012
).
50.
M.
Srebro
and
J.
Autschbach
,
J. Phys. Chem. Lett.
3
,
576
(
2012
).
51.
I. E.
Brumboiu
,
G.
Prokopiou
,
L.
Kronik
, and
B.
Brena
,
J. Chem. Phys.
147
,
044301
(
2017
).
52.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
53.
J.
Hubbard
,
Proc. R. Soc. A
276
,
238
(
1963
).
54.
V. L.
Campo
, Jr.
and
M.
Cococcioni
,
J. Phys.: Condens. Matter
22
,
055602
(
2010
).
55.
J.
Hubbard
,
Proc. R. Soc. A
285
,
542
(
1965
).
56.
J.
Hubbard
,
Proc. R. Soc. A
296
,
82
(
1967
).
57.
H. J.
Kulik
and
N.
Marzari
,
J. Chem. Phys.
134
,
094103
(
2011
).
58.
H. J.
Kulik
and
N.
Marzari
,
J. Chem. Phys.
135
,
194105
(
2011
).
59.
T.
Körzdörfer
,
J. S.
Sears
,
C.
Sutton
, and
J.-L.
Brédas
,
J. Chem. Phys.
135
,
204107
(
2011
).
60.
S. R.
Whittleton
,
X. A.
Sosa Vazquez
,
C. M.
Isborn
, and
E. R.
Johnson
,
J. Chem. Phys.
142
,
184106
(
2015
).
61.
V.
Vlček
,
H. R.
Eisenberg
,
G.
Steinle-Neumann
,
L.
Kronik
, and
R.
Baer
,
J. Chem. Phys.
142
,
034107
(
2015
).
62.
T. Z. H.
Gani
and
H. J.
Kulik
,
J. Chem. Theory Comput.
12
,
5931
(
2016
).
63.
G.
Moynihan
,
G.
Teobaldi
, and
D. D.
O’Regan
,
Phys. Rev. B
94
,
220104
(
2016
).
64.
J. D.
Gledhill
,
F.
De Proft
, and
D. J.
Tozer
,
J. Chem. Theory Comput.
12
,
4879
(
2016
).
65.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A.
Dal Corso
,
S.
de Gironcoli
,
S.
Fabris
,
G.
Fratesi
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitonsen
,
A.
Smogunov
,
P.
Umari
, and
R.
Wentzcovitch
,
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
66.
A. D.
Becke
,
J. Chem. Phys.
138
,
074109
(
2013
).

Supplementary Material