Laplacian-dependent models of the kinetic energy density: Applications in subsystem density functional theory with meta-generalized gradient approximation functionals
aa r X i v : . [ c ond - m a t . o t h e r] F e b Laplacian-dependent models of the kinetic energy density: Applications in subsystemdensity functional theory with meta-generalized gradient approximation functionals
Szymon ´Smiga,
1, 2, 3
Eduardo Fabiano,
4, 3
Lucian A. Constantin, and Fabio Della Sala
4, 3 Istituto Nanoscienze-CNR, Italy Institute of Physics, Faculty of Physics, Astronomy and Informatics,Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland Center for Biomolecular Nanotechnologies @UNILE,Istituto Italiano di Tecnologia (IIT), Via Barsanti, 73010 Arnesano (LE), Italy Institute for Microelectronics and Microsystems (CNR-IMM),Via Monteroni, Campus Unisalento, 73100 Lecce, Italy
The development of semilocal models for the kinetic energy density (KED) is an important topicin density functional theory (DFT). This is especially true for subsystem DFT, where these modelsare necessary to construct the required non-additive embedding contributions. In particular, thesemodels can also be efficiently employed to replace the exact KED in meta-Generalized Gradient Ap-proximation (meta-GGA) exchange-correlation functionals allowing to extend the subsystem DFTapplicability to the meta-GGA level of theory.Here, we present a two-dimensional scan of semilocal KED models as linear functionals of thereduced gradient and of the reduced Laplacian, for atoms and weakly-bound molecular systems. Wefind that several models can perform well but in any case the Laplacian contribution is extremelyimportant to model the local features of the KED. Indeed a simple model constructed as the sum ofThomas-Fermi KED and 1/6 of the Laplacian of the density yields the best accuracy for atoms andweakly-bound molecular systems. These KED models are tested within subsystem DFT with variousmeta-GGA exchange-correlation functionals for non-bonded systems, showing a good accuracy ofthe method.
I. INTRODUCTION
Advances in Kohn-Sham density functional theory[1, 2] are strictly related to the availability of accurateand efficient approximations for the exchange-correlation(XC) functional. Over the years, many XC functionalshave been proposed [3, 4]. The simplest ones, namelythe local density approximation (LDA) [1] and the gen-eralized gradient approximations [5] (GGAs), display anexplicit dependence on the electron density (and its gra-dient). Thus, they are practical tools within the KSscheme. More advanced methods are developed intro-ducing additional input ingredients, that are not explicitfunctionals of the density. The most notable examples inthis context are the hybrid [6, 7] and the meta-GGA [8, 9]XC functionals. The former include a fraction of Hartree-Fock exchange energy and are thus true non-local ap-proaches. The meta-GGA functionals instead have thegeneral form E xc = Z e xc ( ρ ( r ) , ∇ ρ ( r ) , τ KS ( r )) d r , (1)where τ KS ( r ) = 12 occ . X i |∇ φ i ( r ) | (2)is the positive defined kinetic energy density (KED) (with φ i being the occupied KS orbitals). This ingredient,which characterizes the meta-GGA level of theory, al-lows the functionals to distinguish different density re-gions and to achieve high accuracy for a broad range of systems and properties [8–30]. For these reasons meta-GGA XC functionals are gaining increasing popularity[8, 10–12, 14, 15, 17–21, 27, 29–33]. For a recent Reviewof meta-GGAs functionals, see Ref. 9.The positive defined kinetic energy density is a semilo-cal quantity which is much easier to compute and managethan a fully non-local one such as the exact exchange en-ergy, but it is not an explicit functional of the electrondensity. Thus, meta-GGA functionals cannot be straight-forwardly used in the KS scheme [34–36]. Instead, theyare usually implemented in the generalized KS (GKS)framework [37]. This allows to exploit in the best waythe semilocal nature of the functional and obtain a fa-vorable computational cost. However, it may pose somedifficulties since, in contrast to the KS one, the GKS XCpotential is not local [9, 36].One relevant case where the non-local nature of thepotential is particularly troublesome is the subsystemformulation of DFT [38–40], where the nearsightednessprinciple [41] is employed to partition a complex systemsinto simpler and smaller subsystems whose interaction isdescribed by embedding potentials. In fact, within sub-system DFT the total ground-state electron density ρ ofa given system can be written as ρ ( r ) = M X I =1 ρ I ( r ) ; ρ I ( r ) = N I X i =1 (cid:12)(cid:12) φ Ii ( r ) (cid:12)(cid:12) , (3)where M is the number of subsystems, N I is the numberof occupied KS orbitals in the I -th subsystem and φ Ii arethe KS orbitals of the I -th subsystem, which are obtainedby solving the coupled equations [42, 43] (cid:20) − ∇ + v Is ( r ) + v Iemb ( r ) (cid:21) φ Ii ( r ) = ǫ Ii φ Ii ( r ) . (4)Here v Is is the ordinary KS potential relative to subsystem I , whereas v Iemb = X J = I v JH + δT nadds δρ I + δE naddxc δρ I (5)is the embedding potential which takes into accountthe presence of the other subsystems. Equation (5)includes the classical Coulomb potential v H (electron-nuclei and electron-electron terms) of all the other sub-systems. Moreover, it contains the functional derivativesof the non-additive kinetic and XC functionals T nadds = T s " M X I =1 ρ I − M X I =1 T s [ ρ I ] (6) E naddxc = E xc " M X I =1 ρ I − M X I =1 E xc [ ρ I ] . (7)The non-additive terms of Eqs. (6) and (7) as well as therelated functional derivatives cannot be calculated in adirect manner for functionals whose dependence on thedensity is not known explicitly. This is always the casefor kinetic energy, which therefore requires an appropri-ate semilocal approximation [44–51]. No problem arisesinstead for LDA or GGA XC functionals. To go beyondthe GGA level of theory, however, the standard subsys-tem DFT formalism based on the KS scheme cannot beeasily applied.Recently, several advances of subsystem DFT havebeen proposed such as the use of subsystems’ frac-tional occupations [52], the use of orbital-dependent func-tionals [53], and the embedding of wave-functions inDFT [54, 55]. Moreover, particularly interesting for thepresent work are an extension to the GKS framework[56] and applications of subsystem DFT calculations us-ing meta-GGA XC functionals [57, 58]. In Ref. 57, wepresented a formal approach for subsystem DFT withmeta-GGAs: the operational equations are analogous toEq. (4) but v s is replaced by the GKS potential of thesubsystem while the XC embedding contribution [Eq. (7)and the related derivative] is approximated by a propersemilocal expression. This is obtained by substituting,in the meta-GGA XC functional [Eq. (1)], the orbital-dependent KED τ KS with an approximated semilocalmodel τ . This computational approach has proven to beremarkably accurate, in spite of the well known difficultyof describing the kinetic energy with semilocal approx-imations and the fact that the XC functional dependsnon-linearly on the KED. The reason may be tracedback to the fact that in non-additive contributions aprominent role is played by valence regions, whereas themore complicated core contributions to the kinetic energy are not much relevant. Nevertheless, in previous stud-ies only one meta-GGA XC functional was considered[the Tao-Perdew-Staroverov-Scuseria [8] (TPSS) one] to-gether with two simple models for the kinetic energy den-sity. Thus, a deep understanding of the methodology isstill lacking.In this paper we aim at improving this work. There-fore, in Section II we introduce a whole family of semilo-cal KED models and in Section IV we test them on dif-ferent properties and systems, including subsystem DFTcalculations using different meta-GGA XC functionals.In this way, we can obtain not only a full assessment ofthe method and of the various KED models but also amore complete understanding of the basic features thatare required to construct successful semilocal KED mod-els. Note that we restrict our study to quite simple KEDmodels, displaying a simple dependence on the gradientand a linear one on the Laplacian of the density. Thischoice allows to obtain reasonably accurate results avoid-ing, at the same time, an excessive complexity in themodels that would prevent the possibility of a detailedanalysis and of a quite complete understanding of theunderlying physics. Of course, in search of very accurateand realistic KED models one is required, in general, togo beyond simple functional forms as those consideredhere. However, this increases significantly the complex-ity of the problem also because the exploration of thecorresponding huge fitting space calls for the use of ad-vanced computational techniques (e.g. machine learn-ing). Hence, we leave this task for future work. II. SEMILOCAL KED MODELS ANDMETHODOLOGY
Most of semilocal models for τ KS can be written as τ = τ TF ( F ( s ) + b q ) , (8)where τ TF = (3 / π ) / ρ / is the Thomas-Fermi(TF) KED [59–61], s = |∇ ρ | / [2(3 π ) / ρ / ] is the re-duced gradient, q = ∇ ρ/ [4(3 π ) / ρ / ] is the reducedLaplacian, F ( s ) is a GGA enhancement factor (i.e. afunction of s ), and b is a coefficient. Functionals with thegeneral form in Eq. (8) are named Laplacian-Level meta-GGA (LLMGGA) kinetic energy functionals. If b = 0,instead, we have simple GGA kinetic energy functionals.Among the KED approximations based on Eq. (8),we mention some of particular relevance for the presentwork:- The simple TF plus von Weizs¨acker approxima-tion (TFW) obtained setting b = 0 and F ( s ) =1 + F W ( s ) with F W = (5 / s , so that τ W = τ TF F W ( s ) is the von Weizs¨acker KED [62]. TheTFW (GGA) model describes correctly the density-tail asymptotic region in different applications [63].However, it does not yield in general a very accu-rate KED [57].- The second-order gradient expansion (GE2) [64,65], with F ( s ) = F GE2 ( s ) = 1 + µ GE2 s , µ GE2 =5 /
27 and b = 20 /
9. The GE2 (LLMGGA) describescorrectly the slowly-varying density limit.- The Yang’s general formula, obtained from the one-body Green’s function in the mean-path approxi-mation considering the Feynman path-integral ap-proach [66]. This has F ( s ) = F Y ( s ) = 1 + 5 − b s , (9)where b is the coefficient in Eq. (8). The Yang’sLLMGGA model describes accurately the KED ofatoms and molecules [67], in particular when b =10 /
9, which corresponds to F Y ( s ) = F GE2 ( s ).- The modified GE2 (MGE2) [68] approximation,with F ( s ) = F MGE2 ( s ) = 1 + µ MGE2 s , µ MGE2 =0 . b = 20 /
9. The MGE2 (LLMGGA) de-scribes accurately large neutral atoms.- The APBEK and revAPBEK GGA functionals [47,48] with b = 0 and enhancement factor F [ rev ] AP BEK ( s ) = 1 + (cid:18) κ + 1 µ MGE2 s (cid:19) − , (10)where κ = 0 .
804 for APBEK and κ = 1 .
245 for re-vAPBEK. These functionals, in particular revAP-BEK, give very accurate non-additive kinetic em-bedding energies for non-covalently interacting sys-tems [47, 48].- The τ L LLMGGA model, having F ( s ) = F revAPBEk ( s ) and b = 20 /
9. This is an extensionof the revAPBEK functional, designed to yield thesame kinetic energy but a better description of theKED. Indeed, it has been shown to perform well, insubsystem DFT calculations, when it is used to re-place τ KS inside a TPSS meta-GGA XC functional[57].Note that, the Laplacian term in Eq. (8) gives no con-tribution to the total kinetic energy (for finite systems)nor to the kinetic potential; see appendix A. Thus, it isoften omitted in many applications. However, in Refs.57, 67, 69 it has been shown that the Laplacian term is avery important quantity to describe accurately the spa-tial dependence of the KED. Nevertheless, because of thepresence of this term, most of the models based on Eq.(8) violate one or both of the following exact conditions: τ KS ≥ τ KS ≥ τ W . (12)In fact, as q can be negative in some regions of space,e.g. near the nucleus, then Eq. (11) is easily violated if b = 0. Instead, the condition in Eq. (12) is harder tosatisfy and it is actually respected only by few accurate functionals [70–73] and by the simple TFW approxima-tion. To avoid this drawback, models based on Eq. (8)can be regularized forcing τ to recover τ W near the nu-cleus [67, 70] (see also Eqs. (19) and (20) later on). Inthis way, however, the Laplacian term will not integrateto zero, thus the regularization will change the kineticenergy and its functional derivative.To assess the performance of these models and espe-cially to understand the role of the gradient and Lapla-cian contributions therein, in this work we consider thegeneral ansatz τ ( a, b ) = τ T F (cid:0) as + bq (cid:1) , (13)where a and b are parameters. Equation (13) provides aquite flexible expression which recovers most of the afore-mentioned KED models. On the other hand, Eq. (13)does not include high-order terms in s , such as the onein Eq. (10). Nevertheless, typical meta-GGA exchangefunctionals, i.e. TPSS, show a significant dependence onthe KED only for small values of s , see Ref. [8]. Thus,the KED needs to be well approximated only in slowly-varying density region.We perform a full two-dimensional scan of the param-eters a and b in Eq. (13) in order to study the possiblebest performance of the ansatz and understand the roleof the different contributions depending on s and q . Tothis purpose, we consider, for each system, the followingerror indicators (functions of a and b ) δ z ( a, b ) = 1 N Z ρ ( r ) (cid:12)(cid:12) z KS − z ( a, b ) (cid:12)(cid:12) d r , (14) δ α ( a, b ) = 1 N Z ρ ( r ) (cid:12)(cid:12) α KS − α ( a, b ) (cid:12)(cid:12) d r , (15) δ x ( a, b ) = 1 N Z ρ ( r ) (cid:12)(cid:12) F x ( s, z KS , α KS ) − F x ( s, z ( a, b ) , α ( a, b )) (cid:12)(cid:12) d r , (16)where N is the number of electrons of the system, F x denotes the TPSS exchange enhancement factor [8], and z KS = τ W /τ KS , (17) α KS = ( τ KS − τ W ) /τ TF , (18) z ( a, b ) = 1 − max (cid:0) , τ ( a, b ) − τ W (cid:1) τ ( a, b ) , (19) α ( a, b ) = max (cid:0) , τ ( a, b ) − τ W (cid:1) τ TF . (20)The quantities z KS and α KS are the main τ KS -dependentvariables used in the construction of the meta-GGA ex-change [9]. In Eqs. (19) and (20) they are regularized inorder to satisfy the constraint of Eqs. (11) and (12). Aspreviously mentioned, the regularization is very impor-tant near the core if b = 0; see also section Sec. IV A.The indicators defined in Eqs. (14)-(16) express howmuch an approximated KED impacts on average on theaccuracy of z KS and α KS as well as on a typical exchangeenhancement factor. Thus, they provide a measure onthe expectable accuracy of each model for subsystemDFT calculations.Finally, to obtain an overall assessment, we average,for each indicator, over M different systems as∆ θ ( a, b ) = 1 M M X i δ θ ( a, b ) , (21)where θ = z, α, or x . III. COMPUTATIONAL DETAILSA. Two-dimensional scan
To study the model of Eq. (13), we have performed atwo dimensional scan over the parameters space ( a, b ) ∈ [0 , × [0 , ,Ne-Ar,(CH ) (WI); (H S) , CH Cl-HCl, CH SH-HCN(DI); (HF) , (H O) , HF-HCN (HB); AlH-HCl,LiH-HCl,BeH -HF (DHB); NH -F , C H -ClF, H O-ClF (CT).For this set we have calculated the error indicators δ z , δ α , and δ x defined by Eqs. (14)-(16). All quantities havebeen calculated on densities obtained from standard KS-DFT supermolecular calculations, using the PBE [74]XC functional. In some cases, aprroximate KED havebeen regularized imposing the constraint of Eq. (12), i.e. τ approx ( reg. ) = max( τ approx , τ W ); see Eqs. (19), (20),and (31).All molecular calculations have been performed usinga locally modified version of the TURBOMOLE [75, 76]program package. In order to guarantee a good accu-racy of the results and to minimize numerical errors, thesupermolecular def2-TZVPPD [77, 78] basis set was em-ployed in all calculations, together with very accurate in-tegration grids ( grid 7 , radsize 14 ). For all KS-DFTcalculations tight convergence criteria were enforced, cor-responding to maximum deviations in density matrix el-ements of 10 − a.u.For the noble atoms we used a fully numerical code[63] and the LDA functional. B. Subsystem DFT
To assess the performance of the various τ models, wehave employed them to carry on subsystem DFT calcu-lations with meta-GGA XC functionals. In this case, wehave computed the non-additive XC meta-GGA termsusing the TPSS [8], revTPSS [10], BLOC [12, 13, 79],meta-VT { } [80], and MGGA MS2 [20] functionals.Note that the latter functional uses a GGA correlationexpression, so the approximation concerns in this caseonly the exchange term. For all subsystem DFT calcula-tions, we have considered a partition in two subsystemsand we have performed a full relaxation of embeddedground-state electron densities through freeze-and-thaw cycles [38, 43]. These calculations have been performedusing the FDE script [56] of the TURBOMOLE program[75] considering as convergence criterion the difference ofthe dipole moments of the embedded subsystems ( < − a.u.). To compute the non-additive kinetic contributions[Eq. (6)] the revAPBEK kinetic functional [47, 48] hasbeen employed. As to the non-additive meta-GGA XCterm, we have employed the computational procedure de-scribed in Ref. 57 (i.e. substitution of τ KS with a semilo-cal model) using the τ models described in Section II.Such calculations have been performed on the followingsets of non-covalent complexes: WI : weak interaction (He-Ne, He-Ar, Ne , Ne-Ar, CH -Ne, C H -Ne, (CH ) ) DI : dipole-dipole interaction ((H S) , (HCl) , H S-HCl, CH Cl-HCl,CH SH-HCN, CH SH-HCl) HB : hydrogen bond ((NH ) , (HF) , (H O) , HF-HCN,(HCONH ) , (HCOOH) ) DHB : double hydrogen bond (AlH-HCl, AlH-HF, LiH-HCl, LiH-HF, MgH -HCl, MgH -HF, BeH -HCl,BeH -HF) CT : charge transfer (NF -HCN,C H -F ,NF -HCN,C H -Cl , NH -F , NH -ClF, NF -HF, C H -ClF,HCN-ClF, NH -Cl , H O-ClF, NH -ClF).The geometries as well as the reference binding energieswere taken from Refs. 46, 81–84.The accuracy of the subsystem DFT calculations hasbeen measured considering the following two quantities.The first one is the error on the total embedding energy∆ E , defined as∆ E = E A + B [˜ ρ A , ˜ ρ B ] − E GKS [ ρ GKS ] , (22)where E A + B [˜ ρ A , ˜ ρ B ] is the total energy obtained fromthe subsystem DFT calculation, with ˜ ρ A and ˜ ρ B beingthe embedded subsystem densities and E GKS being theconventional GKS total energy of the complex, computedat the “true” ground-state density ρ GKS [48, 56, 85]. Thesecond one is the embedding density error defined as ξ = 1000 N Z | ∆ ρ ( r ) | d r , (23)where N is the total number of electrons and ∆ ρ ( r ) isthe deformation density∆ ρ ( r ) = ˜ ρ A ( r ) + ˜ ρ B ( r ) − ρ GKS ( r ) . (24)Note that only valence electron densities have been con-sidered in the calculation of density errors [57].Finally, to have an overall indication of the perfor-mance of the different τ approximations, we computed,within each group of molecules, the the mean abso-lute error (MAE) and the mean absolute relative error(MARE), the latter being referred to reference bindingenergies [85]. Furthermore, we considered the two globalquantities [48]rwMAE = 15 X i (cid:18) M AE i < M AE i > (cid:19) (25)rwMARE = 15 X i (cid:18) M ARE i < M ARE i > (cid:19) (26)where the sums extend over WI, DI, HB, DHB, CT, and < M AE i > ( < M ARE i > ) is the average MAE (MARE)among the different methods considered for the class ofsystems i . C. Embedding energy error decomposition
The error on embedding energy can be written [57, 85]∆ E = ∆ D + ∆ T S + ∆ E xc . (27)The error components are∆ D = E H "X I ρ I + e E xc "X I ρ I + (28) − ( E H (cid:2) ρ GKS (cid:3) + e E xc "X I ρ GKS , ∆ T S = e T nadds − T S [Φ GKS ] − X I T s [Φ I ] , (29)∆ E xc = e E xc [ ρ GKS ] − X I e E xc [ ρ I ] − ( E xc [Φ GKS ] − X I E xc [Φ I ] ) , (30)where E H is the Hartree energy (electron-electron andelectron-nuclei classical Coulomb interaction energy), ρ GKS and Φ
GKS are respectively the density and theSlater determinant of the supermolecular system as ob-tained from a standard GKS calculation, and the tilde( e ) denotes that an approximate functional form is con-sidered. IV. RESULTSA. Two dimensional scan of the τ model Figure 1 reports the values of the error indicators δ z , δ α , and δ x [see Eqs. (14)-(16)] as computed for the Neonatom in the case of the KED model of Eq. (13). Similarresults are found for Argon and Radon, see Fig. S1 ofRef. [86].The most important result of Fig. 1 is that the minimafor all the three indicators appear at a = 0 and b ≈ .
2: the latter almost coincide with the “non-empirical”value of 20/9 (from the second-order gradient expansion).
FIG. 1: Values of δ z ( a, b ) , δ α ( a, b ) and δ x ( a, b ) [see Eqs. (14)-(16)] as functions of the values of the a and b parametersin Eq. (13) for the Neon atom. The white horizontal linesdenote a = 5 /
27 and a = 5 /
3. The white vertical lines denote b = 10 / b = 20 /
9. The white diagonal line representsthe Yang’s formula [Eq. (9)]. The color map is linear and,in each panel, it is limited to represent values not larger thanfour times the minimum value.
Similar findings have been obtained in Ref. [69], wherethe value of b had been fitted for different GGA KED, inorder to reproduce the exact KED for the first 10 atoms(H-Ne) of the periodic table. The result in Fig. 1, is moregeneral, because we demonstrate that ( a, b ) = (0 , .
2) isthe global minimum.Thus we introduce a new KED model, named regular-ized Thomas-Fermi plus Laplacian (TFL), that is definedby the formulas τ T F L ( reg. ) = max (cid:0) τ T F L , τ W (cid:1) , (31) τ T F L = τ T F (cid:18) q (cid:19) . (32)We note again that the TFL model can be considered a“non-empirical” functional because it is defined with thecoefficient b = 20 / τ W . However, the regularization is only active near thecore (see hereafter), which is not relevant for molecularapplications. Thus, for simplicity, we will consider bothTFL and TFL(reg.) as LLDA functionals, to distinguishthem from other LLMGGA functionals which depend onthe gradient also in the valence region.Fig. 1 shows that the shapes of δ z and δ x are rathersimilar. This similarity is indeed not surprising because,as mentioned before, the z ingredient, tested in δ z , is themain variable used to construct the meta-GGA TPSSexchange enhancement factor, which is used in the con-struction of δ x . On the other hand, the δ α plot looksquite different. Nevertheless, the minimum is still at( a, b ) ≈ (0 , / a, b ) = (0 , a, b ) = (5 / , δ z and δ α . This is due to errorcancellation effects, as explained in Fig. 2, where theintegrands of Eq. (14) and Eq. (15) are reported in theupper and lower panel, respectively.Concerning the integrand of δ z (upper panel), wesee that TFW always underestimates the exact value,whereas TF(reg.) overestimates it, yielding similar globalerror δ z . Concerning the integrand of δ α (lower panel),one can see that TF(reg.) is better than TFW in therange to r < < r < < r < r > δ α . Clearly, without regularization the results forTF will be much worse (see red dashed-lines in Fig. 2). -5051015 z π r ρ [ a . u . ] KSTF (reg.)TFWTFL (reg.)TFTFL r [a.u.] -5051015 α π r ρ [ a . u . ] FIG. 2: Integrand of Eq. (14), upper panel, and Eq. (15),lower panel, versus the radial distance for the Neon atom fordifferent KED models.
Figure 2 also shows the high accuracy of the TFL(reg.)model which matches the exact KS results everywhere inthe space. Note that the regularization is very impor-tant when z KS needs to be accurately reproduced nearthe core region, otherwise, as shown in the upper panelof Fig. 2, a divergence (i.e. at r ≈ .
02) will appear, as q becomes negative near the core. On the other hand,regularization is less critical for α KS , as one can see com-paring TFL and TFL(reg.) in the lower panel Fig. 2.The above results are not limited to atoms, but holdalso for molecules. Figure 3 reports a similar investiga-tion as in Fig. 1 but averaged over a set of molecules(see Section III). Interestingly, the results for the threeindicators are all similar to the ones shown in Fig. 1,with the minima again located at ( a, b ) ≈ (0 , / B. Comparison with other kinetic energyfunctionals
In Table I the performances of KED from different ki-netic energy functionals are compared to TFL: we consid-ered LDA, GGA and LLMGGA functionals. Note thatall functionals are regularized, but TFW, mGGA andmGGArev which implicitly satisfy Eq. (12). The resultsin Tab. I show that TF and all GGAs give almost thesame results for all indicators, as already discussed insection IV A.Large improvements are only obtained when the Lapla-cian term is considered. The mGGA functional of Ref.[70] yields the best results for ∆ z and ∆ x , while TW02 TABLE I: Values of ∆ z ( a, b ) , ∆ α ( a, b ) and ∆ x ( a, b ) [see Eq.(21)] for different KED models evaluated on a molecular test set(see Section III for details). If applicable, also the values of the a and b parameters (as defined in Eq. (13)) and the literaturereference are reported for each functional. The smallest values of the indicators for each class of functionals are highlighted inbold style.Model Reg. ∆ z ∆ α ∆ x a b Ref.Local Density Approximation (LDA)TF Yes 0.193 0.662 0.017 0 0 [61]Generalized Gradient Approximation (GGA)TFW Impl. τ TF (1 + µ GE2 s ) Yes 0.189 τ TF (1 + µ MGE2 s ) Yes 0.188 τ L Yes 0.070 0.241 0.008 - 20/9 [57] τ TW02 + τ TF q Yes 0.068 τ LC94 + τ TF q Yes 0.068 0.236 0.007 - 20/9 [88]mGGA Impl. - - [70]L0.4 Yes 0.189 0.663 0.017 - - [49]L0.6 Yes 0.189 0.663 0.017 - - [49]mGGArev Impl. 0.096 0.323 0.010 - - [73]Laplacian-Level Density Approximation (LLDA)TFL Yes [87] regularized with (20 / q , gives the best results for∆ α ,The last line of Tab. I shows that TFL model yieldsthe best results overall, confirming the results of Sect.IV A. C. Subsystem DFT calculations
In this section we apply the TFW, τ L , and τ T F L [Eq.(32)] KED models in subsystem DFT calculations usingseveral different XC meta-GGA functionals. This canprovide a practical assessment of the reliability of thedifferent models. We consider TFW as a simple GGAsatisfying Eq. (12), τ L as a simple meta-GGA, previ-ously considered in Ref. [57], and τ T F L as the best func-tional, as discussed in the previous sections. Note that,in subsystem DFT, τ L and τ T F L can be also used with-out regularization . In fact, as discussed in Ref. [57], corecontributions cancel out in the non-additive quantities.In Fig. 4 we report (panel a) the rwMAE [Eq. (25)]for the embedding densities errors and (panel b) the rw-MARE [Eq. (26)] for the embedding energy errors eval-uated for all the investigated meta-GGA XC functionalsand KED models. Results for all systems are reported inRef. 86. Concerning the density error, we see that the τ models containing the Laplacian term (i.e. τ L and τ T F L ) perform very similarly and accurately, yieldingrwMAE values in the range 0 . − .
98. The τ T F L isbetter than τ L only for the MS2 functional, while theperformances are comparable for other meta-GGAs. Onthe other hand, the τ T F W model displays a worse per-formance, providing rwMAE values oscillating between1 . − . τ L and τ T F L provide the bestperformance with quite similar results, but for the BLOCfunctional where τ T F L is distinctively better than τ L . Inall the cases the TFW model gives definitely poorer re-sults yielding rwMARE values above 1.1 for all function-als (to be compared to values of about 0.9 for the τ T F L and τ L models) and it is very inaccurate for revTPSSand MS2.For a deeper understanding of these results, it is pos-sible to perform an embedding energy error decomposi-tion analysis [57, 85] (see subsection III C for details).This allows to separate the relaxation error ∆ D , due tothe fact that the embedding density differs from the ref-erence Kohn-Sham supermolecular one, from the errorsarising from the approximations used in the non-additiveenergy functionals. The latter can be further divided intoa kinetic energy error ∆ T S and an XC error ∆ E xc . The aaa bbb δ α δ x δ z FIG. 3: Values of ∆ z ( a, b ) , ∆ α ( a, b ) and ∆ x ( a, b ) [see Eqs.21)] as functions of the values of the a and b parameters inEq. (13) for a set of molecules (see Sect. III). Other detailsare the same as in Fig. 1. TPSS revTPSS BLOC meta-VT{8,4} MS2 r w M A R E TFW τ L τ TFL
TPSS revTPSS BLOC meta-VT{8,4} MS2 r w M A R E TFW τ L τ TFL a)b)
FIG. 4: Top: the rwMAE [Eq. (25)] for the embedding den-sities errors calculated for all the investigated meta-GGA XCfunctionals and τ models. Bottom: the rwMARE [Eq. (26)]for the embedding energy errors calculated for all the investi-gated meta-GGA XC functionals and τ models. former originates from approximations used in Eq. (6)and, in the present calculations, it is the same for all cases(we used the revAPBEK kinetic functional in all calcu-lations). The latter traces back to the use of a semilocalKED model to compute Eq. (7), for different meta-GGAXC functionals. Thus, this is the most interesting quan-tity in the present context.Hence, we consider in Fig. 5 the average XC absoluteratio (AXCAR) [57] defined asAXCAR = 1 M M X i =1 | ∆ E xc || ∆ E xc | + | ∆ T s + ∆ D | , (33)where M is the number of systems in the test set (resultsfor all systems are reported in supporting information).This provides an average absolute (i.e. without errorcompensation) measure of the XC contribution to theembedding energy error. Inspection of the figure showsthat in all cases τ L and τ T F L yield AXCAR values muchsmaller than TFW (AXCAR alwayes larger than 40%).Moreover, τ T F L gives smaller AXCAR than τ L , for allfunctionals, but revTPSS. These results indicate that theaccuracy observed for the embedding energy error is notessentially due to some error compensation effect, butreally traces back to a superior quality of the τ T F L modelin approximating the kinetic energy density contributionsin the XC meta-GGA non-additive term.
V. SUMMARY AND CONCLUSIONS
In subsystem DFT the non-additive XC term must bean explicit functional of the density. Therefore, in or-der to be able to perform subsystem DFT calculations
TPSS revTPSS BLOC meta-VT{8,4} MS2 AX C A R % TFW τ L τ TFL
FIG. 5: The AXCAR error [Eq. (33)] calculated for differentmethods and τ models. employing meta-GGA XC functionals, it is necessary toreplace the positive-defined Kohn-Sham KED ( τ KS ) withan approximated semilocal KED model ( τ ). The devel-opment of such a model is not as a hard task as thegeneral development of accurate kinetic energy function-als because, in the context of subsystem DFT, only non-additive contributions are important. Thus, core con-tributions, which incorporate most of the Pauli kineticenergy, are not much relevant. On the other hand, in thepresent context, there is the need to model all the spatialfeatures of the KED, since this is employed non-linearlyinside meta-GGA XC functionals and not only to pro-duce a kinetic energy (via integration) and/or a kineticpotential (via functional derivation). For this reason, forexample, the accurate determination of the proper gaugefor the KED is a crucial issue in this field. The devel-opment of accurate KED models is therefore a challeng-ing and interesting topic that has practical relevance insubsystem DFT but is also of fundamental importanceto understand the basic behavior of the kinetic energy,which is a key quantity in all areas of DFT.In this work, we studied the problem of developing anaccurate semilocal KED model by considering a flexibleansatz [Eq. (13)] depending linearly on the reduced gra-dient and Laplacian of the density. In this way, we foundthat a Laplacian contribution of the type (1 / ∇ ρ isessential to fix the correct gauge in the KED and thusto provide an accurate description of most of the fea-tures of the exact KDE. In fact, already a simple model(i.e. the TFL one of Eq. (31)) built from the Thomas-Fermi KED and the Laplacian correction (i.e. withoutany gradient correction), performs remarkably well foratoms and molecular complexes.A thorough assessment work showed that the simpleTFL model, is indeed able to capture most of the impor-tant features of the KED. Hence, in this context it out-performs all the known gradient expansion KED models(e.g. the conventional [64], modified [68], path-integral-derived [66] second-order gradient expansions) as well asmany GGA and Laplacian-dependent meta-GGA mod-els.These findings are important for different topics inDFT. Clearly they are relevant for subsystem DFT cal-culations with meta-GGA XC functionals. Thus, theycan allow to extend the applicability and accuracy ofsubsystem DFT, permitting to benefit from the advan-tages of meta-GGA XC functionals. Nevertheless, thegeneral results of this work, concerning the role of dif- ferent density contributions to the KED and especiallythat of the Laplacian, can have a larger impact on fu-ture work. They are in fact extremely relevant for allthose studies and applications where a semilocal modelof the KED is important. Here we suggest, for example,the evaluation of local indicators such as the electron lo-calization function (ELF) [89, 90] or the entanglementlength [91]. These are density indicators but they are ac-tually obtained as orbital-dependent expressions due tothe presence of τ KS . The use of an appropriate semilocalKED model can turn them into true density indicatorsand largely extend their application domain providing aswell a deeper understanding of the underlying physicsand chemistry. Acknowledgements
This work was partially supported by the Na-tional Science Center under Grant No. DEC-2013/11/B/ST4/00771.
Appendix A: Total energy and potential of themodel KED of Eq. (8)
Consider the KED of Eq. (8). Note that it can bewritten τ = τ TF F ( s ) + 3 b ∇ ρ = τ GGA + 3 b ∇ ρ . (A1)The corresponding total kinetic energy is T s = Z τ GGA d r + 3 b Z ∇ ρd r = T GGAs + 3 b Z ∇ ρd r . (A2)The last integral can be evaluated via the second Green’sidentity. Hence, Z ∇ ρd r = Z ∇ ⊥ ρdS + Z ρ ∇ d r = Z ∇ ⊥ ρdS , (A3)where ∇ ⊥ ρ is the perpendicular component of the gradi-ent with respect to the surface area element dS . For anyfinite system, with an exponentially decaying density, theintegral in Eq. (A3) vanishes. Therefore, the Laplaciantem does not contribute to the total kinetic energy.Concerning the potential we have [9] v T = ∂τ∂ρ − ∇ · ∂τ∂ ∇ ρ + ∇ ∂τ∂ ∇ ρ . (A4)For τ given by Eq. (8), that is v T = ∂τ GGA ∂ρ − ∇ · ∂τ GGA ∂ ∇ ρ + ∇ ∂ (3 b ∇ ρ/ ∂ ∇ ρ = (A5)= v GGAT + 3 b ∇ ∂ ∇ ρ∂ ∇ ρ = v GGAT + 3 b ∇ v GGAT . [1] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[2] J. F. Dobson, G. Vignale, and M. P. Das, Electronic Den-sity Functional Theory (Springer, 1998).[3] G. E. Scuseria and V. N. Staroverov, in
Theory andApplications of Computational Chemistry: The First 40Years (A Volume of Technical and Historical Perspec-tives) , edited by C. E. Dykstra, G. Frenking, K. S. Kim,and G. E. Scuseria (Elsevier, Amsterdam, 2005), pp. 669–724.[4] J. P. Perdew and K. Schmidt, AIP Conf. Proc. , 1(2001).[5] D. C. Langreth and M. J. Mehl, Phys. Rev. B , 1809(1983).[6] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem.Phys. , 9982 (1996).[7] A. V. Arbuznikov, J. Struct. Chem. , S1 (2007).[8] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuse-ria, Phys. Rev. Lett. , 146401 (2003).[9] F. Della Sala, E. Fabiano, and L. A. Constantin, Int. J.Quantum Chem. , 1641 (2016).[10] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Con-stantin, and J. Sun, Phys. Rev. Lett. , 026403 (2009).[11] L. A. Constantin, L. Chiodo, E. Fabiano, I. Bodrenko,and F. Della Sala, Phys. Rev. B , 045126 (2011).[12] L. A. Constantin, E. Fabiano, and F. Della Sala, Phys.Rev. B , 035130 (2012).[13] L. A. Constantin, E. Fabiano, and F. Della Sala, Phys.Rev. B , 125112 (2013).[14] F. Della Sala, E. Fabiano, and L. A. Constantin, Phys.Rev. B , 035126 (2015).[15] T. Van Voorhis and G. E. Scuseria, J. Chem. Phys. ,400 (1998).[16] H. L. Schmider and A. D. Becke, J. Chem. Phys. ,8188 (1998).[17] Y. Zhao and D. G. Truhlar, J. Chem. Phys. , 194101(2006).[18] A. Ruzsinszky, J. Sun, B. Xiao, and G. I. Csonka, J.Chem. Theory. Comput. , 2078 (2012).[19] J. Sun, B. Xiao, and A. Ruzsinszky, J. Chem. Phys. ,051101 (2012).[20] J. Sun, R. Haunschild, B. Xiao, I. W. Bulik, G. E. Scuse-ria, and J. P. Perdew, J. Chem. Phys. , 044113 (2013).[21] R. Peverati and D. G. Truhlar, J. Phys. Chem. Letters , 117 (2012).[22] B. Xiao, J. Sun, A. Ruzsinszky, J. Feng, R. Haunschild,G. E. Scuseria, and J. P. Perdew, Phys. Rev. B ,184103 (2013).[23] J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao,A. Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P.Perdew, Phys. Rev. Lett. , 106401 (2013).[24] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P.Perdew, Phys. Rev. B , 075102 (2004).[25] C. Adamo, M. Ernzerhof, and G. E. Scuseria, J. Chem.Phys. , 2643 (2000).[26] K. E. Riley, B. T. Op’t Holt, and K. M. Merz, J. Chem.Theory. Comput. , 407 (2007).[27] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. , 036402 (2015).[28] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky,H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, et al.,Nat. Chem. , 831 (2016), article. [29] H. S. Yu, X. He, and D. G. Truhlar, J. Chem. Theory.Comput. , 1280 (2016).[30] J. Tao and Y. Mo, Phys. Rev. Lett. , 073001 (2016).[31] J. Wellendorff, K. T. Lundgaard, K. W. Jacobsen, andT. Bligaard, J. Chem. Phys. , 144107 (2014).[32] N. Mardirossian and M. Head-Gordon, J. Chem. Phys. , 074111 (pages 1) (2015).[33] L. A. Constantin, E. Fabiano, J. Pitarke, andF. Della Sala, Phys. Rev. B , 115127 (2016).[34] S. K¨ummel and L. Kronik, Rev. Mod. Phys. , 3 (2008).[35] A. V. Arbuznikov and M. Kaupp, Chem. Phys. Lett. ,495 (2003).[36] F. Zahariev, S. S. Leang, and M. S. Gordon, J. Chem.Phys. , 244108 (2013).[37] A. Seidl, A. G¨orling, P. Vogl, J. A. Majewski, andM. Levy, Phys. Rev. B , 3764 (1996).[38] T. A. Wesolowski, in Chemistry: Reviews of CurrentTrends , edited by J. Leszczynski (World Scientific: Sin-gapore, 2006, Singapore, 2006), vol. 10, pp. 1–82.[39] C. R. Jacob and J. Neugebauer, Wiley Interdisci-plinary Reviews: Computational Molecular Science ,325 (2014).[40] A. Krishtal, D. Sinha, A. Genova, and M. Pavanello, J.Phys.-Condens. Mat. , 183202 (2015).[41] E. Prodan and W. Kohn, PNAS , 11635 (2005).[42] T. A. Wesolowski and A. Warshel, J. Phys. Chem. ,8050 (1993).[43] T. A. Wesolowski and J. Weber, Chem. Phys. Lett. ,71 (1996).[44] A. W. G¨otz, S. M. Beyhan, and L. Visscher, J. Chem.Theory Comput. , 3161 (2009).[45] T. A. Wesolowski, J. Chem. Phys. , 8516 (1997).[46] T. A. Wesolowski, H. Chermette, and J. Weber, J. Chem.Phys. , 9182 (1996).[47] L. A. Constantin, E. Fabiano, S. Laricchia, and F. DellaSala, Phys. Rev. Lett. , 186406 (2011).[48] S. Laricchia, E. Fabiano, L. A. Constantin, and F. DellaSala, J. Chem. Theory. Comput. , 2439 (2011).[49] S. Laricchia, L. A. Constantin, E. Fabiano, and F. DellaSala, J. Chem. Theory. Comput. , 164 (2014).[50] F. Tran and T. A. Wesoowski, Int. J. Quantum Chem. , 441 (2002).[51] F. Tran and T. A. Wesolowski, in Recent Progress inOrbital-free Density Functional Theory , edited by T. A.Wesolowsky and Y. A. Wang (World Scientific, Singa-pore, 2013), pp. 429–442.[52] E. Fabiano, S. Laricchia, and F. Della Sala, J. Chem.Phys. , 114101 (2014).[53] S. Laricchia, E. Fabiano, and F. Della Sala, Chem. Phys.Lett. , 114 (2011).[54] T. A. Weso lowski, Phys. Rev. A , 012504 (2008).[55] T. Dresselhaus and J. Neugebauer, Theor. Chem. Acc. , 1 (2015).[56] S. Laricchia, E. Fabiano, and F. Della Sala, J. Chem.Phys. , 164111 (2010).[57] S. ´Smiga, E. Fabiano, S. Laricchia, L. A. Constantin, andF. Della Sala, J. Chem. Phys. , 154121 (2015).[58] P. Ramos, M. Papadakis, and M. Pavanello, J. Phys.Chem. B , 7541 (2015).[59] L. H. Thomas, Proc. Cambridge Phil. Soc. , 542(1926). [60] E. Fermi, Rend. Accad. Naz. Lincei , 73 (1928).[61] E. Fermi, Z. Phys. , 602 (1927).[62] C. F. von Weizs¨acker, Z. Phys. A , 431 (1935).[63] C. Cirac`ı and F. Della Sala, Phys. Rev. B , 205405(2016).[64] D. A. Kirzhnitz, Field theoretical methods in many-bodysystems (Pergamon Press, 1967).[65] M. Brack, B. Jennings, and Y. Chu, Phys. Lett. B , 1(1976).[66] W. Yang, Phys. Rev. A , 4575 (1986).[67] W. Yang, R. G. Parr, and C. Lee, Phys. Rev. A , 4586(1986).[68] D. Lee, L. A. Constantin, J. P. Perdew, and K. Burke, J.Chem. Phys. , 034107 (2009).[69] D. Garc`ıa-Aldea and J. E. Alvarellos, J. Chem. Phys. , 144109 (2007).[70] J. P. Perdew and L. A. Constantin, Phys. Rev. B ,155109 (2007).[71] V. V. Karasiev, R. S. Jones, S. B. Trickey, and F. E.Harris, Phys. Rev. B , 245120 (2009).[72] V. V. Karasiev, D. Chakraborty, O. A. Shukruto, andS. B. Trickey, Phys. Rev. B , 161108 (2013).[73] A. C. Cancio, D. Stewart, and A. Kuna, J. Chem. Phys. , 084107 (2016).[74] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[75] TURBOMOLE V6.2, 2009, a development of Universityof Karlsruhe and Forschungszentrum Karlsruhe GmbH,1989-2007, TURBOMOLE GmbH, since 2007; availablefrom .[76] F. Furche, R. Ahlrichs, C. H¨attig, W. Klopper, M. Sierka,and F. Weigend, Wiley Interdisciplinary Reviews: Com- putational Molecular Science , 91 (2014).[77] F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. , 3297 (2005).[78] D. Rappoport and F. Furche, J. Chem. Phys. ,134105 (2010).[79] L. A. Constantin, E. Fabiano, and F. Della Sala, J. Chem.Theory. Comput. , 2256 (2013).[80] J. M. del Campo, J. L. G´azquez, S. Trickey, and A. Vela,Chem. Phys. Lett. , 179 (2012).[81] Y. Zhao and D. G. Truhlar, J. Phys. Chem. A , 5656(2005).[82] Y. Zhao and D. G. Truhlar, J. Chem. Theory and Com-put. , 415 (2005).[83] S. Laricchia, E. Fabiano, and F. Della Sala, J. Chem.Phys. , 124112 (2013).[84] E. Fabiano, L. A. Constantin, and F. Della Sala, J. Chem.Theory. Comput. , 3151 (2014).[85] S. Laricchia, E. Fabiano, and F. Della Sala, J. Chem.Phys. , 014102 (2012).[86] See supplementary material athttp://aip.scitation.org/doi/suppl/10.1063/1.4975092.[87] F. Tran and T. A. Wesoloski, Int. J. Quant. Chem. ,441 (2002).[88] A. Lembarki and H. Chermette, Phys. Rev. A , 5328(1994).[89] A. D. Becke and K. E. Edgecombe, J. Chem. Phys. ,5397 (1990).[90] B. Silvi and A. Savin, Nature , 683 (1994).[91] S. Pittalis, F. Troiani, C. A. Rozzi, and G. Vignale, Phys.Rev. B91