The Role of the Reduced Laplacian Renormalization in the Kinetic Energy Functional Development

: The Laplacian of the electronic density diverges at the nuclear cusp, which complicates the development of Laplacian-level meta-GGA (LLMGGA) kinetic energy functionals for all-electron calculations. Here, we investigate some Laplacian renormalization methods, which avoid this divergence. We developed two different LLMGGA functionals, which improve the kinetic energy or the kinetic potential. We test these KE functionals in the context of Frozen-Density-Embedding (FDE), for a large palette of non-covalently interacting molecular systems. These functionals improve over the present state-of-the-art LLMGGA functionals for the FDE calculations.


Introduction
The non-interacting kinetic energy (KE) is a central quantity in different electronic structure theories, such as orbital-free density functional theory (OF-DFT) [1][2][3], subsystem density-functional theory (DFT) and frozen-density-embedding [4,5], and hydrodynamic models [6][7][8]. Despite the fact that the development of accurate KE functional is an old problem, it has received only very small interest, if compared, e.g., to the one devoted for exchange-correlation (XC) functionals [9,10].
More recently, the class of Laplacian-level meta-GGA (LLMGGA) KE functionals, in which the Laplacian of the density (∇ 2 ρ) is used as an additional ingredient, has attracted strong interest [11,[23][24][25][26][27][28][29]. The incorporation of the latter is especially important because it allows for distinguishing different density regions (e.g., the nuclear and the bonding region), which are out of reach for all present day GGAs. Furthermore, with this ingredient, we are able to satisfy additional constraints [23,30], such as the recovery of the fourth-order gradient expansion (GE4) [23,25,31,32], or its modifications [33], leading to more robust KE functionals. Laplacian-level meta-GGA functionals have also been investigated for XC functionals [23,[34][35][36][37][38][39] showing that the use of ∇ 2 ρ in the functional development, may lead to some difficulties, such as the appearance of unphysical oscillations in the corresponding potential near the nucleus. Clearly, this problem is not present when pseudopotentials are used [11,27].
The near-core divergence was investigated by Cancio et al. [39], in the context of XC functional construction. It was suggested that a renormalization of the Laplacian term may lead to a significant improvement in the quality of XC potentials. Therefore, in the current work, we will investigate this idea in the context of KE functional development. We also propose a simpler and more accurate renormalization of the Laplacian.
In order to study the methods for the improvement of the quality of the meta-GGA kinetic potential, we will build very simple KE functionals, with enhancement factors similar to the modAPBE exchange functional one [39], by using, to the full extent, the conjointness conjecture hypothesis [40,41]. We construct two similar LLMGGA functionals that practically differ only from the methods used to obtain a smooth kinetic potential (i.e., renormalization procedure of the Laplacian). We show that these functionals are of interest for the Fronzen Density Embedding (FDE) calculations [42][43][44][45][46], significantly improving the FDE performance of meta-GGA KE functionals [25].
The Laplacian renormalization procedure proposed in this work can be further applied not only for the development of KE functionals designed especially for orbital-free DFT [11], but also to the Laplacian-level meta-GGA XC functionals that are of recent interest [45][46][47] because they restore the meta-GGA XC potential to the simpler Kohn-Sham DFT, instead of using the generalized Kohn-Sham scheme.
This article is organized as follows. In Section 2, we construct the modAPBEz and modAPBEq LLMGGA KE functionals, using two strategies for the Laplacian renormalization. The modAPBEz and modAPBEq enhancement factors and their functional derivatives are also discussed in detail. In Section 3, we present the computational details of the FDE and other molecular and atomic calculations that are the subject of Section 4. Finally, in Section 5, we summarize our results and provide conclusions and perspectives.

Theory
The semi-local Laplacian-level KE functional can be written in the general form as where is the approximated kinetic energy density (KED) defined as product of Thomas-Fermi (TF) KED [48][49][50] (τ TF = 3 10 (3π 2 ) 2/3 ρ 5/3 ) and the KE enhancement factor F s , expressed through dimensionless reduced gradient and Laplacian of density respectively. Both ingredients contain information allowing for distinguishing different density regions [51]. This is shown in Figure 1, where we present the spatial behavior of both quantities for two representative systems, namely the Ne atom and the Ne dimer. In the density-tail (where s and q → ∞) and slowly varying regions, both quantities are rather proportional to each other. The qualitatively different behavior is more pronounced near the nuclei and in the bond region. In the latter (see the right panel in Figure 1), the reduced gradient goes to zero, whereas q = 0.
Similar differences can be noted in the nuclear region. The reduced Laplacian tends here to q → −∞, while the reduced gradient takes a small, finite value (s ≈ 0.4). The divergence of q may be especially useful to model the deenhancement effect of F s near the nucleus, which is actually out of reach for standard GGA functionals accurate for the uniform electron gas scaling [52], yielding in this region F s ≥ 1. We recall that, in large atomic systems, the enhancement factor near nuclei tends to [30] F s (s = 0.3534, q → −∞) = 0.2367 . (4) This feature can be easily recovered by imposing certain restrictions on the enhancement factor form [39].  The description of the approximated KED is not sufficient to model an accurate KE functional. In fact, for applications in OF-DFT and subsystem DFT, the most important quantity to be described is the kinetic potential, which, in the case of the LLMGGA level-of-theory, has the form v k (r) = ∂τ app ∂ρ − ∇ · ∂τ app ∂∇ρ The first term in Equation (6) is the TF kinetic potential multiplied by the enhancement factor itself. The second term describes the small correction coming from the variation of F s w.r.t. the density. The third term is the GGA term, and it can cause some difficulties in GGA functional development i.e., divergences at the nucleus [53] and an unphysical behavior of the corresponding potential in the middle of the bond [54]. Furthermore, it imposes certain conditions on the enhancement factor for which 1 s everywhere in the space. This condition is not satisfied by some KE functionals, e.g., Reference [18]. Finally, the last term in Equation (6) comes from a variation with respect to ∇ 2 ρ. This part of the kinetic potential gives rise to major problems, such as unphysical oscillations near the nucleus, or some numerical instability problems related to the high-degree derivatives involved in it (up to ∇ 4 ρ).
As was shown by Cancio et al. [39], the enhancement factor features at the nuclear region can be modeled successfully using a "hybridization" variable z(s, q) = as 2 + bq .
This inhomogeneity parameter naturally appears in a gradient expansion of the exchange hole, in the regime of a slowly varying density [39,55]. It has been used in the context of XC [56] (with a = 40 and b = 10 9 ) and KE [26] (with a = − 40 27 and b = 20 9 ) functional development to recover the correct behavior of functionals around the nuclei. Here, the a and b parameters in Equation (8) were set to 8 30 and b = 0.2, according to Reference [39]. When the new "hybridization" variable z is used, the condition in Equation (7) can be rewritten as 1 s ∂F s ∂z ∂z ∂s = 2a ∂F s ∂z < ∞. In the following, we present two new KE functionals, which use a renormalized reduced Laplacian.

ModAPBEz
In order to construct a KE functional using Equation (8), we start with the modPBE [39] enhancement factor F modPBE where µ 0 , η and κ are parameters. Equation (9) guarantees the finite behavior of the enhancement factor when z diverges. However, it may produce strong unphysical oscillations in the corresponding potential v k (r) = δT s [ρ]/δρ(r) (see the discussion in Section 2.4), due to the ∇ 2 ∂τ app ∂∇ 2 ρ term. As it was pointed out in Reference [39], this problem can be partially solved by utilization of a "renormalized" version of z [26,39] where C describes the decay rate of the exponential function and H(x) is a Heaviside function. Both z and z are presented in Figure 1. One can note that, while z exhibits the same features as q in the vicinity of nuclei, the z, due to renormalization, approaches the finite negative value −Cκ/(3µ 0 ). Following Reference [57], we also introduce the "renormalized" µ 0 as which gives additional control of the kinetic potential curvature and converges to µ 0 (1 − A) for z → −∞. The final enhancement factor, named modAPBEz, is: which has the following limits The first limit allows for fixing the µ 0 parameter, in order to recover the modified second-order gradient expansion (MGE2) value of µ MGE2 = 0.23889 [15]. This can be easily done by substituting Equation (8) into this limit, whereby we find µ 0 = µ MGE2 /0.8. The last limit can be combined with the Equation (4), to fix the enhancement factor at the nuclei and thus we obtain the relation The remaining parameters C, η, and κ can be chosen through optimization of kinetic potential curvature, in the nuclear region, where we minimize the following integral [39] quantity for few atomic systems. The optimization of κ (which is usually responsible for the tail asymptotic behavior of the enhancement factor) together with C can be justified by its strong correlation with η value (κ(η)). It was shown in Reference [39] that the variation of both parameters influences the asymptotic behavior of the enhancement factor given by Equation (9). Thus, setting η = 3 [39], we perform only the optimization of the remaining two parameters. The optimization of Equation (14) was performed with respect to the He, Ne, and Kr densities obtained from analytical Hartree-Fock orbitals [58]. For each atom, the value of Equation (14) was evaluated (I k ). As a threshold for optimization, we have considered two conditions: These conditions have been met for C = 0.26839, and κ = 4.0147 so that A = 0.634054 . The κ parameter is approximately five times larger than κ 0 used in the APBEκ [15] energy functional. Hence, the construction of the modAPBEz meta-GGA KE functional is completed.

ModAPBEq
Instead of performing a renormalization of the whole z (as in Equation (10)), we can renormalize only the part that actually causes the divergence at the nuclei, namely the reduced Laplacian ingredient. For this purpose, we define a "renormalized" q as with the following properties: We recall that a quite similar ingredient has been used in the construction of the Laplacian-level EDMGGA exchange functional [56]. In Figure 1, we plotted the q r for the Ne atom and the Ne dimer. In the tail of density as well as in the middle of the bond, q and q r behave qualitatively the same.
Then, we define a new "renormalized" z as z qr (s, q r ) = as 2 + bq r , with a and b kept unchanged. One can note, looking at Figure 1, that z and z qr behave almost identical in the whole region of space, being superimposed. The major difference only appears in the nuclear region where both variables tend to different finite values, namely z min = −1.20438 and z min qr = −0.0666956. Now, substituting z → z qr in Equation (9), we obtain the new KE functional (denoted here as modAPBEq), which has the following limits 1+3ηµ 0 z min qr /κ+(3µ 0 z min qr /κ) 2 , z qr → z min qr .
As previously, µ 0 was set to recover µ in the slowly varying density limit, and we fixed again η = 3. Thus, only κ is left as a free parameter, which we finally fixed to κ = 3.216 optimizing w.r.t. subsystem DFT calculations.
Even if modAPBEq does not recover Equation (4), its behavior at the nucleus is still accurate, as shown in the next subsection. Finally, we remark that both modAPBEz and modAPBEq meta-GGA functionals recover the MGE2 in the slowly-varying density limit. This fact insures that they behave correctly under the uniform-electron-gas density scaling [52], and under the Thomas-Fermi density scaling [52,59], then being accurate in the semi-classical limit of large, neural, and non-relativistic atoms [60].
For simplicity, henceforth we will use for modAPBEz and modAPBEq the shortened names mAPBEz and mAPBEq, respectively.

Comparison of KE Enhancement Factors
In Figure 2, we present enhancement factors as functions of the reduced Laplacian q (left panels) and reduced gradient s (right panels) for several LLMGGA functionals. The plots are reported for two values of q (q = 0 and q = 2) and s (nuclei (s = 0.3534) and bond (s = 0)). The density regime characterized by s = 0.3534 and negative q is important in the nuclear region. The L0.4 [25] LLMGGA KE functional shows a large enhancement in this region, which is incorrect from the physical point of view. On the other hand, the mAPBEz goes to the exact deenhancement of Equation (4) and mAPBEq is also less than one.
A similar behavior of enhancement factor can be observed in the bond region where (s = 0) and q can be either positive or negative (depending on the bond type [61,62]). All GGA KE functionals give here F s = 1. The L0.4 exhibits a rather strong enhancement for any type of bonding. In case of mAPBEz and mAPBEq functionals, we can note deenhancement for the bond types characterized by the negative value of q.
In the case where q = 0 (right panels), and s ≤ 1, all functionals behave similarly, recovering either MGE2 (as modAPBEz and modAPBEq) or GE4 (as L0.4). For larger values of s (s > 1), the functionals start to vary from each other, reaching (for large values of s) the finite value determined by the κ parameter.
Finally, for a positive value of q (q = 2), the shape of L0.4 and mAPBEq enhancement factors for s < 0.6 are rather similar, while the mAPBEz functional shows a stronger enhancement.

Comparison of KE Potentials
As a final element of this section, we discuss the impact of the inclusion of the Laplacian of the density on the quality of the kinetic potentials.
In order to check the overall behavior of the potentials generated by mAPBEz and mAPBEq functionals, we report them in Figure 3 for two representative systems, namely the He and Ne atoms, respectively. Additionally, for comparison, we also plotted the MGE2 and L0.4 [25] potentials. All potentials have been calculated using analytical Hartree-Fock orbitals [58]. The exact potentials are taken from Reference [2]. The figures clearly show a rather poor behavior of all potentials in comparison to "exact" reference values. Clearly, the errors are larger for the smaller systems, the He atom, for which the exact KE functional is just the von Weizsäcker [63] one, which is not considered in the construction of the present meta-GGA functionals, which have instead increased accuracy with the dimensions of the systems (Recall that MGE2 is exact for an atom with an infinite number of electrons [15]).
The errors are clearly visible both in the core and in long-range asymptotic limit, where the exact kinetic potential approaches the energy of the highest occupied molecular orbitals, in contrast to all meta-GGA functionals considered here which decay to zero. We recall that, in order to improve the overall quality of the kinetic potential, the functional must have a positive Pauli enhancement factor and potential [11,64]. However, by construction of mAPBEz, mAPBEq, and L0.4 meta-GGAs functionals have PBE-like enhancement factors, and, consequently, the von Weizsäcker [63] expression is not recovered asymptotically. However, we recall that our main interest in the present study is the smoothness of the potentials near the core.
In Figure 3, we observe that the mAPBEq potentials are much smoother than the L0.4 and mAPBEz ones, which have large oscillations close to the nucleus (in the case of the He atom, the L0.4 and mAPBEz go to negative values). Instead, the mAPBEq potential is remarkably accurate and as smooth as the MGE2 potential. This fact clearly shows that the renormalized ingredient q r (see Equation (15)) has important features, being much simpler but also superior to the whole renormalization and smoothness procedures used in the construction of mAPBEz.

Computational Details
To enable the comparison of the quality of new functionals, we performed the following benchmark calculations: •

Subsystem DFT calculations:
We considered a partition of the total density into two fragments ρ = ρ A + ρ B , where ρ A and ρ B are densities corresponding to subsystems A and B, respectively. The full relaxation of embedded ground-state electron densities was obtained using the freeze-and-thaw cycles [42,43] and considering convergence when the difference of dipole moments of the embedded subsystems is below 10 −3 a.u. In case of KS-DFT calculation, the maximum deviations in density matrix elements of 10 −7 a.u. were considered as convergence criteria. The benchmark set consists of five weakly interacting groups of molecular complexes used in our previous studies [21,45,46].
We considered the embedding density error (ξ) and the total embedding energy error (∆E). The first is defined as where N is the total number of electrons and ∆ρ(r) = ρ A (r) + ρ B (r) − ρ KS (r) is the deformation density calculated as the difference between the converged subsystem densities (ρ A , ρ B ) and the conventional ground state KS density (ρ KS ). In Equation (18), only the valence electron density was considered.
The ∆E is defined as where E A+B [ρ A , ρ B ] is the total energy obtained from the subsystem DFT calculation, with ρ A and ρ B being the embedded subsystem densities, and E KS is the conventional KS total energy of the complex, computed for ground-state density ρ KS [15,65].
All calculations have been performed using the PBE [70] exchange-correlation functional, computed with the TURBOMOLE program package [71].

Subsystem DFT Calculations
The newly defined functionals were employed in subsystem DFT calculations to compute the non-additive KE part of the total energy. The performance of the KE functionals was evaluated using the mean absolute error (MAE) that was calculated for each group of complexes separately. Moreover, to investigate overall performance of KE functionals, we have calculated in the same spirit the relative weighted MAE (rwMAE) [15] defined as where the sum runs over all group of complexes and < MAE i > is the average MAE calculated for all functionals within each group of systems i. In order to compare our data, we have conducted calculations with several state-of-the-art GGA functionals: revAPBEκ [15,16], LGAP [21] and one Laplacian-level functional, namely L0.4 [25]. In case of other Laplacian-level KE functionals such as MGE4 or MGGA [23], the use of a monomolecular subsystem DFT approach is required, in order to guarantee good convergence. As was pointed out in Reference [25], the supermolecular approach gives rise to oscillations in the kinetic potential in the tail of the density, therefore causing convergence problems. Hence, those results are not reported in present work. For the same reasons, we do not report in this section the results obtained with the functional proposed in Reference [26], which exhibits the same problems.
In Table 1, we report the ξ errors. First of all, we note that all calculations conducted with mAPBEz and mAPBEq KE functionals provide stable FDE solutions. This is not the case for L0.4 meta-GGA functional where the lack of convergence is observed for five CT complexes. This is the consequence of the aforementioned oscillations [25] in the kinetic potential. Despite this fact, the MAE and rwMAE suggest that the overall performance of both GGAs and meta-GGA KE functionals is rather similar (the rwMAE range between 0.91-1.14). The relatively similar ξ errors suggest that the inclusion of the Laplacian ingredient does not outperform the GGA functionals. This is probably related to the significant cancellation of errors in the non-additive meta-GGA KE functional and potential [25], which are quite similar to those obtained for regular GGAs (see the discussion of Section 2.4). However, at the Laplacian-level of theory, the newly proposed modAPBEq and modAPBEz functionals are better than the L0.4 functional.
Let us turn our attention to the ∆E errors that are reported in Table 2. Note that the best overall results are obtained using the mAPBEq functional, which yields a rwMAE of 0.90. It outperforms (in three groups of complexes, namely DI, DHB, and CT) the state-of-the-art revAPBEk and LGAP GGA functionals that both give a total rwMAE of 0.94. The mAPBEq results are closely followed by the ones of the mAPBEz functional, and both overcome the L0.4 functional. Interestingly, in case of DHB complexes, the best MAE (1.14) is provided by the L0.4 KE functional being almost twice as small as those reported by other functionals, but, for all other cases, L0.4 has lower accuracy. In conclusion, mAPBEz and mAPBEq improve over the current state-of-the-art Laplacian-level L0.4 functionals leading to a stable convergence in FDE calculation and providing the results that are in line with current state-of-the-art GGAs.

Results for Atoms and Molecules
In Table 3, we show the results for small atoms and molecules. All functionals give similar errors for the IP13 and EA13 tests. In this case, the best functional is mAPBEq. For the other tests (aKE, mKE, and atKE), the best results are found from the APBEk GGA, followed by revAPBEk GGA. All meta-GGAs worsen the total KE of atoms and molecules. Even if mAPBEz has been optimized for He, Ne, and Kr atoms, it is worse than the non-empirical GE4 and L0.4, for both aKE and mKE tests. However, for the atomization KE test (atKE), the mAPBEz and mAPBEq are significantly better than GE4 and L0.4, being as accurate as the revAPBEk GGA. This is the most important result of this subsection, showing that mAPBEz and mAPBEq improve over L0.4 not only for FDE calculations, but also for molecular atomization KE.

Conclusions
We have proposed the mAPBEz and mAPBEq Laplacian-level meta-GGA KE functionals, constructed from the mAPBE exchange functional [39], by using the conjointness conjecture hypothesis [40,41]. Both mAPBEz and mAPBEq functionals are similar, and they mainly differ from the renormalization and smoothness procedures used to dump/cancel the non-physical oscillations introduced by the Laplacian of the density, in the kinetic potential close to the nuclear cusp.
However, while mAPBEz uses to the full extent the renormalization and smoothness methods proposed in Reference [39], the mAPBEq is much simpler, being based only on the renormalized reduced Laplacian q r of Equation (15). Moreover, as shown in Figure 3, the mAPBEq kinetic potential is everywhere remarkably smooth. On the other hand, the mAPBEz kinetic potential still has some unwanted structure at the nuclear cusp. These results fully assess the q r ingredient that can and should be of interest for further, Laplacian-dependent KE development.
We have tested the mAPBEz and mAPBEq meta-GGA KE functionals in the context of FDE, for a large palette of non-covalently interacting molecular systems, spanning several important classes of forces, such as weak-interaction, dipole-dipole, hydrogen-bond, double-hydrogen-bond, and charge transfer. These functionals improve and extend the applicability over the L0.4 meta-GGA [25], which was the present state-of-the-art Laplacian-level functional for FDE calculations [25]. The overall accuracy of the mAPBEz and mAPBEq is in line with, or even better than, the revAPBEk GGA [15] and LGAP GGA [21]. Finally, we note that the development of LLMGGA is still in its infancy. GGA functionals have been more deeply investigated, but their higher accuracy is based on a strong error compensation on different density regions. Such error compensation can be largely eliminated using the Laplacian ingredient [25], but this task is more difficult to achieve and different paths need to be investigated in the future.