Theory and Functionals in the Hypercomplex Kohn-Sham DFT Method
TTheory and Functionals in the Hypercomplex Kohn-Sham DFT Method
Neil Qiang Su ∗ Department of Chemistry, Key Laboratory of Advanced Energy Materials Chemistry (Ministry of Education) andRenewable Energy Conversion and Storage Center (RECAST), Nankai University, Tianjin 300071, China
The static/strong correlation error inherent in existing functionals has severely limited the furtherapplication of Kohn-Sham (KS) density functional theory (DFT). The recently developed hyper-complex KS (HCKS) theory shows great potential to overcome the fundamental limitations of theconventional KS-DFT, hence further development and application of HCKS will effectively guidethe construction of new functionals toward better description of strong correlation. To this end,this work derives the working equations for the HCKS calculation and proves that HCKS usingcommon functionals has much lower computational complexity than reduced density matrix func-tional theory (RDMFT), while maintaining the same density search space as RDMFT. Moreover,this work combines HCKS with hybrid functionals that are based on orbitals and their occupations.The test on triplet-singlet gaps shows that the systematic error of HCKS due to the lack of staticcorrelation in approximate functionals can be effectively reduced by including appropriate amountof static correlation correction, while the large systematic errors of RKS and UKS rarely changewith the correction. Therefore, HCKS creates new channels for the development and evaluation ofapproximate functionals, and the insights and advances gained in this work will be helpful to reducethe static correlation error in approximate functionals.
Kohn-Sham (KS) density functional theory (DFT) [1–5] is now the leading electronic structure method for bothmolecular and extended systems. Exact in principle, KS-DFT has to rely on density functional approximations(DFAs) to the exchange-correlation (XC) energy in prac-tice. While existing DFAs have made the great successof KS-DFT, they cause notorious failures in practical ap-plications [6–11], which severely limits the broader appli-cability of DFT.The static/strong correlation issue is one of the mostimportant and long-standing topic in DFT [12–16]. It hasbeen proved that KS-DFT with commonly used DFAs,including local density approximations (LDAs), general-ized gradient approximations (GGAs), meta-GGAs, andhybrid functionals, shows serious static correlation error[5, 10, 11, 17]. Especially for systems where strong cor-relation becomes dominant, i.e. strongly correlated sys-tems such as molecules with stretched bonds and Mottinsulators, KS-DFT with common DFAs utterly fails[10]. Many efforts have been devoted toward system-atic removal of the static correlation error inherent in ap-proximate functionals [17–33]. Work from Becke showedsome promise of describing strong correlation with XChole modeling [30, 34, 35], based on which Johnson andContreras-Garc´ıa constructed strong-correlation modelsfor better description of atoms with fractional chargesand spins [28, 29]. It is worth mentioning that an ef-fective fractional-spin correction was developed recently[17], which can restore the flat-plane behavior of elec-tronic energy and properly describe the chemical bonddissociation without breaking spin symmetry. Neverthe-less, the strong correlation problems that can be solvedare very limited, and the progress to date is still far frombeing able to systematically eliminate the static correla-tion error in approximate functionals. A recent effort proved that the generalization of KS or-bitals to hypercomplex number systems [36] can furtherextend the search space for the density in KS-DFT to aspace that is equivalent to natural spin orbitals with frac-tional occupations in reduced density matrix functionaltheory (RDMFT[37, 38]) . The resulting HypercomplexKS (HCKS)[36] is able to capture the multi-reference na-ture of strong correlation by optimized orbitals with dy-namically varying fractional occupations. HCKS showsgreat potential to overcome the fundamental limitationsof the conventional KS method, which thus creates newchannels for the development and evaluation of approxi-mation functionals.The focus of this work is twofold. On the one hand,The theoretical framework of HCKS will be introducedin detail and extended to general functionals that de-pend on orbitals and their occupation numbers. onthe other hand, this work seeks to gain in-depth in-sights into the development of approximate functionalsthrough the comparison between KS and HCKS, in use ofthe same functionals that combine semi-local functionalswith occupation-dependent Hartree-Fock (HF) like ex-change. The advances in this work should be helpful forboth the development and applications of DFT.The total electronic energy in DFT is uniquely deter-mined by the electronic density E tot [ ρ σ ], which is furtherdivided into the following four terms in the KS theory[2, 3] E tot [ ρ σ ] = T s [ ρ σ ] + E ext [ ρ σ ] + E H [ ρ σ ] + E XC [ ρ σ ] , (1)where E ext [ ρ σ ] is the external energy E ext [ ρ σ ] = α,β (cid:88) σ (cid:90) v ext ( r ) ρ σ ( r ) d r , (2) a r X i v : . [ phy s i c s . c h e m - ph ] F e b E H [ ρ σ ] is the Coulomb energy E H [ ρ σ ] = 12 α,β (cid:88) σσ (cid:48) (cid:90) ρ σ ( r ) ρ σ (cid:48) ( r ) r d r d r , (3)and E XC [ ρ σ ] is the exchange-correlation (XC) energy,which is obtained by approximation in practice. T s [ ρ σ ] isthe kinetic energy of the noninteracting reference system.Because the ground-state wave function of the noninter-acting reference system is just a determinant, i.e. the KSdeterminant, T s [ ρ σ ] can be explicitly expressed in termsof the occupied KS orbitals [2, 3], T s [ ρ σ ] = − α,β (cid:88) σ N σ (cid:88) k =1 (cid:104) ϕ σk |∇ | ϕ σk (cid:105) , (4)while the Kohn-Sham ansatz guarantees that the exactground state density can be written as the ground statedensity of the noninteracting reference, ρ σ ( r ) = N σ (cid:88) k =1 | ϕ σk ( r ) | . (5)The ground-state energy in KS can be obtained by theminimization of E tot [ ρ σ ] with respect to { ϕ σp } , or equiv-alently, by solving the following KS equationsˆ h σ ϕ σp ( r ) = ε σp ϕ σp ( r ) , (6)where the KS operator is defined byˆ h σ = − ∇ + v H ( r ) + v σ XC ( r ) + v ext ( r ) , (7)and v H ( r ) and v σ XC ( r ) are the Coulomb and XC poten-tials, which are obtained via the functional derivatives δE H [ ρ σ ] δρ σ ( r ) and δE XC [ ρ σ ] δρ σ ( r ) respectively.In HCKS theory, the orbitals are formulated in hyper-complex number systems [36], which take the followingform ϕ σp ( r ) = φ σ, p ( r ) + n (cid:88) µ =1 φ σ,µp ( r ) e µ , (8)where { e , e , · · · , e n } are a basis of dimension n in aClifford algebra, such that [39, 40] e µ = − e µ e ν = − e ν e µ , (9)and { φ σ,µp } are a set of real functions. The conjugatehypercomplex of Eq 8 takes the form of¯ ϕ σp ( r ) = φ σ, p ( r ) − n (cid:88) µ =1 φ σ,µp ( r ) e µ . (10)The HCKS orbitals are a set of general KS orbitals of highdimension, and the determinant constructed by occupiedHCKS orbitals is called the HCKS determinant. Although the HCKS orbitals and HCKS determinantare hypercomplex, the resulting density and kinetic en-ergy are real, which can be obtained by inserting Eq 8into Eqs 5 and 4 respectively. They are ρ σ ( r ) = N σ (cid:88) k =1 n (cid:88) µ =0 [ φ σ,µk ( r )] , (11)and T s [ ρ σ ] = − α,β (cid:88) σ N σ (cid:88) k =1 n (cid:88) µ =0 (cid:104) φ σ,µk |∇ | φ σ,µk (cid:105) . (12)Thereby, based on the HCKS orbitals defined by Eq 8,the energy of Eq 1 remains real. The minimization of E tot [ ρ σ ] with respect to the HCKS orbitals leads to aset of high-dimension KS equations, named the HCKSequations, ˆ h σ φ σ, p ( r ) = ε σp φ σ, p ( r )ˆ h σ φ σ, p ( r ) = ε σp φ σ, p ( r ) · · · ˆ h σ φ σ,np ( r ) = ε σp φ σ,np ( r ) . (13)It can be deduced from Eq 13 that all nonzero φ σ,µp for µ = 0 – n are the eigenvectors of ˆ h σ with the same eigen-value ε σp . When the eigenvector associated with ε σp isnot degenerate, the components φ σ,µp in the p -th HCKSorbital are either 0 or different only by a constant fac-tor; while for the degenerate case, degenerate eigenvec-tors and their mix can appear in different componentsof the p -th HCKS orbital, which can lead to fractionaloccupations for each degenerate eigenvectors.Instead of solving Eq 13, { φ σ,µp } can be further ex-panded on a set of orthonormal functions { χ p } φ σ,µp ( r ) = K (cid:88) q =1 χ q ( r ) V σ,µpq , (14)where V σ,µ is a K × K matrix associated with the µ -thcomponent of the HCKS orbitals. With this, the densityof Eq 11 becomes, ρ σ ( r ) = K (cid:88) p,q =1 χ p ( r ) D σpq χ q ( r ) , (15)and D σ is the density matrix on { χ p } D σ = n (cid:88) µ =0 V σ,µT I N σ K V σ,µ , (16)where I N σ K is a K × K diagonal matrix, with the first N σ diagonal elements being 1 and the rest being 0. BecauseD σ is symmetric, it can be diagonalized by an orthogonalmatrix U σ , D σ = U σ Λ σ U σT , (17)where Λ σ is a diagonal matrix, diag( λ σ , λ σ , · · · , λ σK ), and { λ σp } are the eigenvalues of D σ . By inserting Eq 17, thedensity of Eq 15 can be further formulated as ρ σ ( r ) = K (cid:88) p =1 λ σp | χ σp ( r ) | . (18)Here { χ σk } are obtained via the unitary transformationon { χ k } : χ σp ( r ) = (cid:80) Kq =1 χ q ( r ) U σqp , which guarantees that { χ σk } are orthonormal. The kinetic energy in HCKS isobtained by inserting Eqs 14, 16 and 17 into Eq 12, whichtakes the form T s [ ρ σ ] = − α,β (cid:88) σ K (cid:88) p =1 λ σp (cid:104) χ σp |∇ | χ σp (cid:105) . (19)The ground-state energy can be obtained by the mini-mization of E tot [ ρ σ ] with respect to both { χ σp } and { λ σp } ,subject to (cid:104) ψ σp | ψ σq (cid:105) = δ pq ; 0 ≤ n σp ≤ , K (cid:88) p =1 n σp = N σ . (20)In practice, the minimization is performed repeatedly on { χ σp } and { λ σp } respectively, until the changes for both { χ σp } and { λ σp } are small between two iterations. Of bothminimization procedures, the minimization process for { χ σp } is more computational complexity. Different fromRDMFT using costly gradient algorithms, the minimiza-tion with respect to { χ σp } in HCKS can be achieved viathe following eigenvalue equationsˆ h σ χ σp ( r ) = ε σp χ σp ( r ) , (21)which greatly reduces the computational cost of HCKSas compared to RDMFT. Note that although Eq 21 hasthe same form as Eq 6 in use of the same XC functional,HCKS allows dynamically varying fractional occupationswhich can result in a different density from KS, leadingto different ˆ h σ and eigenvectors.The XC functionals in HCKS can be semi-local puredensity functionals, or functionals that implicitly dependon the density through { χ σp } and { λ σp } . Based on { χ σp } and { λ σp } , the Hartree-Fock like (HF-like) exchange be-comes E HFX = − α,β (cid:88) σ N σ (cid:88) k,l n σk n σl (cid:90) χ σk ( r ) χ σ ∗ k ( r ) χ σl ( r ) χ σ ∗ l ( r ) r d r d r . (22)Similarly, the long-range and short-range HF-like ex-change functionals can be defined by including the errorfunction erf( ωr ) and the complementary error functionerfc( ωr ) in Eq 22. With these, the global hybrid func-tionals [41–46] and all kinds of long-range corrected hy-brid functionals [47–60] in KS can be applied in HCKS as well. The ground-state energy of HCKS with hybridfunctionals can be obtained via the similar procedure de-scribed above, except that the XC potential in ˆ h σ of Eq21 becomes delocalized, which includes the HF-like ex-change potential v HF ,σ X ( r , r ) = − N σ (cid:88) k n σk χ σk ( r ) χ σ ∗ k ( r ) r . (23)Similarly, the long-range and short-range parts of the HF-like exchange potential can be defined correspondingly.Here the spin-restricted KS (RKS), spin-unrestricted(UKS), and HCKS in use of the same semi-local and hy-brid functionals will be evaluated on the recent devel-oped TS12 benchmark set [61]. This dataset includes 12data points for triplet-singlet gaps of several atoms anddiatomic molecules. The ground states of these systemsare triplets, while the lowest singlet states are of biradicalcharacter. More tested details about TS12 can be foundin refs 61, 62. The XC functionals applied here are thenonempirical meta-GGA TPSS [63] and the hybrid func-tionals based on TPSS. All calculations were performedusing a local modified version of the NWChem package[64].The test results of the semi-local functional TPSS canbe found in Tab I. As can be seen that TPSS@RKSsystematically overestimates the triplet-singlet gaps ofthe 12 systems in TS12. Due to the limitation of RKS,TPSS@RKS can only provide a closed-shell solution forany singlet state, with all orbitals either doubly occu-pied or empty. Thereby, it cannot capture the multi-reference nature of the singlet biradicals, resulting inoverestimation of the energies of the lowest singlet states.Note that even though RKS can correctly predict thezero spin-densities ( ρ α − ρ β ) for these singlets, it breaksthe space symmetry of the total densities [36, 62]. Un-like TPSS@RKS, TPSS@UKS can predict the biradicalcharacter by breaking the spin symmetry, with two un-paired electrons of opposite spins confined in differentphysical region. However, the destruction of spin sym-metry leads to the over-relaxation of the orbitals occu-pied by the unpaired opposite electrons, thus the en-ergies of the singlet states as well as the triplet-singletgaps are systematically underestimated. Similar resultsand discussions about RKS and UKS can be found inref 62. In contrast, HCKS allows dynamically varyingfractional occupations to capture the multi-reference na-ture of strong correlation. Hence, TPSS@HCKS signif-icantly improves the predicted triplet-singlet gaps, withthe error of TPSS@RKS reduced by half. Even so, thetriplet-singlet gaps are still systematically overestimateddue to the lack of static correlation in TPSS, both MSEand MAE are 7.27 kcal/mol. Actually, this error can bereduced to 2.41 kcal/mol when the PBE XC functional[65] is used. Therefore, approximate functionals to theXC energy has a large influence on the calculation of TABLE I. Triplet-singlet gaps ∆ E T − S (= E S − E T ) (kcal/mol)of various atoms and diatoms from the TS12 benchmark set[61]. Experimental data are used as reference. The XCfunctionals tested here are TPSS and the hybrid functional E hTPSSXC = c HFX E HFX +(1 − c HFX ) E TPSSX + E TPSSC with c X = − − HCKS, which will be further explored below.In order to effectively guide the development of func-tionals toward better description of strong correlation,it is essential to further compare and evaluate the im-pact of different choices of XC functional on RKS, UKSand HCKS when calculating systems of multi-referencenature. To this end, the triplet-singlet gaps of TS12 arecalculated by RKS, UKS, and HCKS in use of the hybridTPSS (named hTPSS) functionals with varying amountof HF exchange ( c HFX ), i.e. E hTPSSXC = c HFX E HFX + (1 − c HFX ) E TPSSX + E TPSSC . It can be clearly seen from Fig1b that hTPSS@RKS systematically overestimates thetriplet-singlet gaps, while hTPSS@UKS underestimatesthem. The systematic overestimation by hTPSS@HCKSat c HFX = 0 (i.e. TPSS@HCKS) becomes serious as c HFX increases, until it reaches the same result as hTPSS@RKSwhen c HFX is large enough. In contrast, the systematicerror of hTPSS@HCKS reduces as c HFX decreases from0, and the smallest MAE is obtained at around c HFX = − c HFX is rewritten in the from: E hTPSSXC = E TPSSXC + ∆ E X , where ∆ E X = − c HFX ( E TPSSX − E HFX ). Itshows that this functional includes a correction ∆ E X tothe TPSS functional. It has been proved that commonlyused semi-local exchange functionals include not only theexchange energy, but also the static correlation effect[68]. Thereby ∆ E X can be regarded as a correction ofstatic correlation to TPSS. By adjusting c HFX to includeappropriate amount of ∆ E X , the systematical error ofTPSS@HCKS can be effectively reduced. Tab I showsthat MSE and MAE of hTPSS@HCKS with c HFX = −
10 5 0 5 10 (a) MAE10 5 0 5 10c
HFX (%)2010010 (b) MSE
RKSUKSHCKS
FIG. 1. Errors of triplet-singlet gaps ∆ E T − S (= E S − E T ) byRKS, UKS and HCKS in use of the hybrid TPSS function-als with different amount of HF exchange ( c HFX ). (a) Meanabsolute errors (MAEs) and (b) mean sign errors (MSEs)are plotted. The tested hybrid functionals take the from: E hTPSSXC = c HFX E HFX + (1 − c HFX ) E TPSSX + E TPSSC . The test setis the TS12 benchmark set [61], and the basis set used is aug-cc-pVQZ [66, 67]. All energies are in kcal/mol. tion ∆ E X has little impact on both RKS and UKS. Thesystematical errors of hTPSS@RKS and hTPSS@UKSrarely change for different c HFX . In fact, similar resultscan be obtained by other hybrid functionals based onLDAs, GGAs, and meta-GGAs. This thus indicates thatthe description of strong correlation poses a great chal-lenge to DFT within the frameworks of RKS and UKS,and explains the slow progress in systematically elimi-nating the strong correlation error inherent in commonlyused DFAs.In summary, this work presented the main equationsfor the calculation of the HCKS theory. The derivationshows that the computational complexity of HCKS canbe greatly reduced by solving a set of eigenvalue equa-tions instead of applying the gradient algorithms that areused in RDMFT, for minimizing the energy with respectto orbitals. Therefore, in used of common functionals,HCKS can have the same search space for the densityas RMDFT, while maintaining the same computationalscaling as KS. In addition, this work further combinedHCKS with hybrid functionals that are based on orbitalsand their occupations. The test on the triplet-singletgaps of the TS12 benchmark set shows that althoughHCKS allows dynamically varying fractional occupationsto capture the multi-reference nature of strong correla-tion, its test results largely depend on the choice of ap-proximate functionals. In particular, the systematic errorof HCKS due to the lack of static correlation in approxi-mate functionals can be effectively reduced by includingappropriate amount of the static correlation correction(i.e. E DFAX − E HFX ). Unlike HCKS, RKS and UKS showeven large systematic errors on the triplet-singlet gaps,worse still, the errors rarely change with the correctionof static correction. Therefore, HCKS provides a feasi-ble approach toward the strong correlation problem, theinsights and advances gained in this work will be help-ful in the development and evaluation of approximationfunctionals.Support from the National Natural Science Foundationof China (Grant 22073049), the Natural Science Founda-tion of Tianjin City (20JCQNJC01760), and Fundamen-tal Research Funds for the Central Universities (NankaiUniversity: No. 63206008) is appreciated. ∗ [email protected][1] P. Hohenberg and W. Kohn, Phys. Rev. , B864(1964).[2] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[3] R. G. Parr and W. Yang, Density-Functional Theoryof Atoms and Molecules , Oxford University Press: NewYork, 1989.[4] R. Dreizler and E. Gross,
Density Functional Theory: AnApproach to the Quantum Many-Body Problem , SpringerBerlin Heidelberg, 2012.[5] N. Mardirossian and M. Head-Gordon, Mol. Phys. ,2315 (2017).[6] A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov,and G. E. Scuseria, J. Chem. Phys. , 194112 (2006).[7] A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov,and G. E. Scuseria, J. Chem. Phys. , 104102 (2007).[8] O. A. Vydrov, G. E. Scuseria, and J. P. Perdew, J. Chem.Phys. , 154109 (2007).[9] P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.Lett. , 146401 (2008).[10] P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.Lett. , 066403 (2009).[11] A. J. Cohen, P. Mori-S´anchez, and W. Yang, Science , 792 (2008).[12] K. Burke, J. Chem. Phys. , 150901 (2012).[13] A. J. Cohen, P. Mori-S´anchez, and W. Yang, Chem. Rev. , 289 (2012).[14] A. D. Becke, J. Chem. Phys. , 18A301 (2014).[15] N. Q. Su and X. Xu, Annu. Rev. Phys. Chem. , 155(2017).[16] J. P. Perdew, A. Ruzsinszky, J. Sun, N. K. Nepal, andA. D. Kaplan, Proc. Natl. Acad. Sci. USA (2021).[17] N. Q. Su, C. Li, and W. Yang, Proc. Natl. Acad. Sci.USA , 9678 (2018).[18] A. J. Cohen, P. Mori-S´anchez, and W. Yang, J. Chem.Phys. , 121104 (2008).[19] P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.A , 042507 (2012).[20] H. van Aggelen, Y. Yang, and W. Yang, Phys. Rev. A , 030501 (2013).[21] S. N. Steinmann and W. Yang, J. Chem. Phys. ,074107 (2013).[22] B. Mussard and J. Toulouse, Mol. Phys. , 161 (2017).[23] R. Haunschild, T. M. Henderson, C. A. Jim´enez-Hoyos,and G. E. Scuseria, J. Chem. Phys. , 134116 (2010). [24] J. J. Phillips, A. A. Kananenka, and D. Zgid, J. Chem.Phys. , 194108 (2015).[25] J. Nafziger and A. Wasserman, J. Chem. Phys. ,234105 (2015).[26] M. Filatov, Wiley Interdiscip. Rev. Comput. Mol. Sci. ,146 (2015).[27] J. Kong and E. Proynov, J. Chem. Theory Comput. ,133 (2016).[28] E. R. Johnson and J. Contreras-Garc´ıa, J. Chem. Phys. , 081103 (2011).[29] E. R. Johnson, J. Chem. Phys. , 074110 (2013).[30] A. D. Becke, J. Chem. Phys. , 074109 (2013).[31] A. Bajaj, J. P. Janet, and H. J. Kulik, J. Chem. Phys. , 191101 (2017).[32] F. Zhou and V. Ozolins, ArXiv e-prints (2017),1710.08973.[33] I. Y. Zhang and X. Xu, Wiley Interdiscip. Rev. Comput.Mol. Sci. , e1490 (2021).[34] A. D. Becke, J. Chem. Phys. , 2972 (2003).[35] A. D. Becke, J. Chem. Phys. , 064101 (2005).[36] N. Q. Su, Unity of kohn-sham density functional the-ory and reduced density matrix functional theory, 2021,arXiv:2102.00394.[37] T. L. Gilbert, Phys. Rev. B , 2111 (1975).[38] M. Levy, Proc. Natl. Acad. Sci. USA , 6062 (1979).[39] I. L. Kantor and A. S. Solodovnikov, Hypercomplex num-bers , Springer-Verlag: Berlin, 1989.[40] J. Vaz Jr. and R. da Rocha Jr.,
An Introduction to Clif-ford Algebras and Spinors , Oxford University Press: NewYork, 2019.[41] A. D. Becke, J. Chem. Phys. , 5648 (1993).[42] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J.Frisch, J. Phys. Chem. , 11623 (1994).[43] M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. ,5029 (1999).[44] C. Adamo and V. Barone, J. Chem. Phys. , 6158(1999).[45] X. Xu and W. A. Goddard, J. Chem. Phys. , 4068(2004).[46] Y. Zhao and D. G. Truhlar, Acc. Chem. Res. , 157(2008).[47] A. Savin, Recent Developments and Applications of Mod-ern Density Functional Theory , Elsevier, Amsterdam,1996.[48] P. M. W. Gill, R. D. Adamson, and J. A. Pople, Mol.Phys. , 1005 (1996).[49] H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem.Phys. , 3540 (2001).[50] T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys.Lett. , 51 (2004).[51] W. Ai, W.-H. Fang, and N. Q. Su, J. Phys. Chem. Lett. , 1207 (2021).[52] O. A. Vydrov, J. Heyd, A. V. Krukau, and G. E. Scuseria,J. Chem. Phys. , 074106 (2006).[53] J. D. Chai and M. Head-Gordon, J. Chem. Phys. ,084106 (2008).[54] J. Heyd and G. E. Scuseria, J. Chem. Phys. , 7274(2004).[55] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem.Phys. , 8207 (2003).[56] R. Baer, E. Livshits, and U. Salzner, Annu. Rev. Phys.Chem. , 85 (2010).[57] T. Stein, H. Eisenberg, L. Kronik, and R. Baer, Phys.Rev. Lett. , 266802 (2010). [58] T. Stein, L. Kronik, and R. Baer, J. Chem. Phys. ,244119 (2009).[59] R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett. ,2810 (2011).[60] L. N. Anderson, M. B. Oviedo, and B. M. Wong, J.Chem. Theory Comput. , 1656 (2017).[61] J. Lee and M. Head-Gordon, J. Chem. Phys. , 244106(2019).[62] J. Lee, L. W. Bertels, D. W. Small, and M. Head-Gordon,Phys. Rev. Lett. , 113001 (2019). [63] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuse-ria, Phys. Rev. Lett. , 146401 (2003).[64] E. Apr´a et al., J. Chem. Phys. , 184102 (2020).[65] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[66] T. H. Dunning, J. Chem. Phys. , 1007 (1989).[67] R. A. Kendall, T. H. Dunning, and R. J. Harrison, J.Chem. Phys. , 6796 (1992).[68] N. C. Handy and A. J. Cohen, Mol. Phys.99