Quasi-dimensional models applied to kinetic and exchange energy density functionals

. We investigate the behavior of three-dimensional 3D exchange energy functional of density-functional theory in anisotropic systems with two-dimensional 2D character and 1D character. The local density approximation (LDA), the generalized gradient approximation (GGA), and the meta-GGA behave as functions of quantum well width. We use the inﬁnite-barrier model (IBM) for the quantum well. In the ﬁrst section, we describe the problem of three-dimensional exchange functional, in the second section we introduce the quasi-2D IBM system, in the third section we introduce the quasi-1D IBM system. Using that an exact-exchange functional provides the correct approach to the true two-dimensional limit, we want to show that the 2D limit can be considered as a constraint on approximate functionals. For the 1D limit case we also propose a new functional obtained with methods completely similar to those of 2D limit.


Introduction
The Kohn-Sham (KS) Density Functional Theory (DFT) [1] is the most used method for electronic structure calculations in quantum chemistry and material science [2][3][4]. In the KSDFT self-consistent scheme, the noninteracting kinetic energy functional T s [n] = dr τ (r) = dr occ j |∇φ j (r)| 2 /2 is treated exactly using the one-particle occupied KS orbitals {φ j (r)}, and only the exchange-correlation (XC) energy functional E xc [n] = dr n(r) xc (r) must be approximated [5]. Here n(r) is the electronic density, xc is the XC energy per particle, and τ is the positive defined kinetic energy density.
Using two simple models, the quasi-2D electron gas and the quasi-1D electron gas we show a fundamental limitation of the 3 nonempirical rungs of the Jacob's ladder, namely the local density approximation LDA and its semilocal extensions, generalized gradient approximation GGA and meta-GGA MGGA, a e-mail: vittoria.urso@gmail.com (corresponding author) the most widely used forms of which are worse than the LDA in the strong 2D limit. An exact-exchange functional provides the correct approach to the true two-dimensional limit. We investigate the performance of three-dimensional density functional E xc [n] in the quasi-two-dimensional electron gas, showing how all three semi-local approximations behave as functions of quantum well width.

Quasi-2D IBM system
The quasi-2D IBM system is defined by the KS potential Then the KS orbitals are plane waves in the (yz)-plane, having the following expression with A and k are the normalization area and the wave vector in the (yz)-plane, and l = 1, 2, 3, .. is the subband index. For the quasi-2D regime, we consider only the lowest level is occupied. Then the density is where r 2D s is the bulk parameter of the 2D uniform electron gas (UEG), which will be kept fixed when L is changed from the maximum value L max [38] to zero. Solving the Schrödinger equation, we find the energy levels When only the lowest level is occupied (E 1,k 2D F < E 2,0 ), the density of states of this system begins to resemble the density of states of a 2D electron gas [39]. Here k 2D F = (2n 2D ) 1/2 is the two-dimensional Fermi wavevector, so Note that r 2D is the 2D bulk Fermi wave vector.

Kinetic energy of the quasi-2DEG
The kinetic energy density (defined by T s = drτ ) is where τ W and τ P are the von Weizsäcker and Pauli kinetic energy densities [40,41]. The first of eq. (7) follows from the following equation where τ T F = (3/10)(k 2D F ) 2 n is the Thomas-Fermi kinetic energy (KE) density, with k 2D F = (3π 2 n) 1 The averaged kinetic energies per particle (defined by T s /N ) are where t 2D s = τ 2D /n 2D is the 2D UEG kinetic energy per particle. Then, the Pauli KE per particle fully recovers the 2D UEG, while the von Weizsäcker part diverges as ∼ L −2 , representing the short-wavelengths oscillations in the x-direction. Noting that and using the Euler equation [42] δT s [n] δn(r) and Eq. (1), we conclude that In Fig. 1, we show T P s /N computed from several KE functionals, versus L/L max . We consider the recently proposed PG1 GGA [41], PGint GGA [43], LKT GGA [44], as well as the popular Thomas-Fermi-Weizsäcker (TFW), second-order gradient expansion (GE2) [45], E00 GGA [46], and the Perdew-Constantin (PC) Laplacian-level meta-GGA [47]. We recall that PG1 and LKT are very accurate for the orbital-free DFT (OFDFT) calculations of bulk solids. On the other hand, PGint is based on the second-order gradient singularity expansion which can mimic the singularity of the jellium linear response at the wave vector k = 2k F . Finally, the PC meta-GGA is a very good model for the kinetic energy density τ , but its functional derivative shows strong unphysical oscillations [48] A reparametrization of the PC KE functional has been proposed in Ref. [49].
The best performance in the quasi-2D regime, is found from PGint GGA, closely followed by PG1 and LKT functionals. The worst performances are given by GE2, E00 GGA, and PC meta-GGA, all of them failing badly in the strong quasi-2D regions (when L/L max → 0). Moreover, GE2 and E00 give wrongly negative T P s /N when L/L max ≤ 0.3. In the upper panel of Fig. 2, we show the Pauli potential v P s = δT P s /δn computed from several KE functionals, using the exact density of the quasi-2D IBM with r 2D s = 2 and L = L max /2. Due to the Euler equation, the exact curve must be a constant representing the kinetic potential of the 2D UEG. Any tested KE functional cannot give a constant Pauli potential. However, LKT and PG1 Pauli potentials have less structure in the region 0.2 ≤ x/L ≤ 0.8. Noting that this is a very hard test for any functional, we consider that LKT and PG1 performances are quite remarkable.
In the lower panel of Fig. 2, we show the averaged Pauli potentialv P s = dx nv P s /N versus L/L max for the case r 2D s = 2. The trend is similar with the one of Fig. 1, with PGint providing the best performance.

Exchange energy of the quasi-2DEG
The first-order density matrix of the quasi-2D IBM is where, without any loss of generality, we chose r = (x, 0, 0), and r = (x , ρ , α) in cylindrical coordinates [50]. Note that n(x) = n 1 (r, r). We also recall that the density matrix of the 2D UEG is n 2D−UEG We calculate the exact exchange energy (EXX) per Fig. 3, we show a comparison between EXX, MVS meta-GGA [51], SCAN meta-GGA [52], and Q2D GGA [53] exchange energies per particle ( x = E x /N ), in the whole quasi-2D regime (0 ≤ L/L max ≤ 1). We recall that all tested semilocal functionals (MVS, SCAN, and Q2D) have been constructed from the quasi-2D condition F x ∝ s −1/2 when s → ∞, that gives a finite value of x = E x /N in the 2D limit [54]. Here F x is the exchange enhancement factor, defined by E The best accuracy is obtained with Q2D GGA that, by construction, recovers the 2D LDA exchange energy per particle, when L → 0. We note that EXX behaves as 2D LDA exchange in the 2D limit, similar with the Pauli kinetic energy that also recovers the 2D LDA kinetic energy. This fact shows that the shortwavelengths oscillations in the x-direction which are essential for the divergence of the von Weizsäcker KE per particle, do not contribute to the exchange energy.

Quasi-1D IBM system
The quasi-1D IBM system is defined by the KS potential where ρ = x 2 + y 2 is the radial distance from the z-axis. The Kohn-Sham orbitals have the form where L z is a normalization length, and φ l (ρ) satisfies the equation whose solutions are of the form φ l (ρ) = AJ 0 ( x 0l L ρ), where x 0l is the lth zero of the Bessel function J 0 . The total energy is The quasi-1D regime is defined by the condition where k 1D F is the 1D Fermi wave vector, which we will keep fixed when we shrink L → 0. Then the quasi-1D IBM KS orbitals are Using the rule In the limit L → 0, the system approaches the 1D UEG.
Varying L is equivalent to performing a non-uniform density scaling in two dimensions, n xy λ (x, y, z) = λ 2 n(λx, λy, z).

R E T R A C T E D A R T I C L E
L max /5, and L = L max /10. Note that N = dρ 2πρn(ρ) = 2k 1D F /π is dependent only on the 1D Fermi wave vector. In the lower panel of Fig. 4, we show the reduced gradient s and Laplacian q for these quasi-1D IBM systems. When ρ → 0 then s ∼ −L −4/3 (k 1D F ) −1/3 ρ+O(ρ 2 ) and q ∼ −0.365/[k 1D F L] 2/3 + O(ρ 2 ). On the other hand, when ρ → L, both s and q are diverging. Finally, we observe that q is almost constant for a large part of the quasi-1D region, and after that increases very sharply.

Kinetic energy of the quasi-1DEG
The quasi-1D IBM kinetic energy density is The averaged kinetic energies per particle are where t 1D s = τ 1D /n 1D is the 1D UEG kinetic energy per particle. Then the Pauli KE per particle fully recovers the 1D UEG, while the von Weizsäcker part diverges as ∼ L −2 , representing the short-wavelengths oscillations in the circular ρ-direction. We observe strong similarities with the quasi-2D case. The von Weizsäcker functional derivative is a constant and using the Euler equation and Eq. (14), we conclude that In Fig. 5, we show the quasi-1D T P s /N versus L/L max , for the same KE functionals used in Fig. 1, when k 1D F = 0.5. We observe that all functionals fail badly, diverging in the 1D limit. The best functional in the strong quasi-1D regime is the PGint GGA, while PC meta-GGA and E00 GGA are relatively accurate in the moderate quasi-1D regime. However, we mention that the quasi-1D IBM is one of the most difficult tests for KE functionals.
In the upper panel of Fig. 6, we show the Pauli potential v P s = δT P s /δn computed from several KE functionals, using the exact density of the quasi-1D IBM with k 1D F = 0.5 and L = L max /2. Due to the Euler equation, the exact curve must be a constant representing the kinetic potential of the 1D UEG. Any tested KE functional can not give a constant Pauli potential. However, LKT and PG1 Pauli potentials have less structure in the region 0 ≤ ρ/L ≤ 0.6. All semilocal KE functionals tested in Fig. 6, have v P → 0 when ρ/L → 1. This feature is a consequence of their behavior in the tail of the density. Finally, we observe close similarities of this figure with the upper panel of Fig. 2.
In the lower panel of Fig. 6, we show the averaged Pauli potentialv P s = dx nv P s /N versus L/L max for the case k 1D F = 0.5. All functionals perform similarly, in the weak and moderate quasi-1D regime (a.i. when L/L max ≥ 0.4), they givev P s quite close to the exact 1D UEG potential, but in the strong quasi-1D regime they fail badly, withv P s diverging when 1D limit is approaching. The trend is similar with the one of Fig. 5, with PGint providing the best performance.

Exchange energy of the quasi-1DEG
The first-order density matrix of the quasi-1D IBM is  n 1 (r, r). On the other hand, the density matrix of the 1D UEG is The exchange energy is

R E T R
where, without any loss of generality, we chose r = (ρ, 0, 0), and we consider the changes of variables t = ρ/L, t = ρ /L, and y = z /L. We note that the 1D UEG exchange energy per particle diverges, because of the known Coulomb divergence in 1D.
Using the non-uniform density scaling in two dimensions, see Eq. (21), we find that a GGA exchange enhancement factor (defined by E x = dr n LDA x F x ) must behave in the quasi-1D regime as This is very different from the quasi-2D case, where F x ∝ s −1/2 , see Ref. [54]. Then we propose the exchange functional, named Q1D GGA, with the following exchange enhancement factor: where a = 0.06525 has been fitted to the quasi-1D IBM. By construction, Q1D GGA is accurate in the quasi-1D density regime, and recovers PBEsol GGA exchange functional at small reduced gradients. In Fig. 7, we show a comparison of several exchange enhancement factors. The Q1D GGA recovers PBEsol only at very small reduced gradients (s ≤ 0.5), and after that F Q1D x sharply decays as a s 2 . In Fig. 8, we report a comparison between EXX, MVS meta-GGA [51], SCAN meta-GGA [52], Q2D GGA [53], and Q1D GGA exchange energies per particle ( x = E x /N ), in the whole quasi-1D regime (0 ≤ L/L max ≤ 1). MVS and SCAN perform almost identical, and Q2D GGA is just a little better, all of them failing badly when L → 0. On the other hand, Q1D GGA is remarkably accurate for the quasi-1D IBM, in both k 1D F = 2 and k 1D F = 0.5 cases. In fact, the same accuracy is obtained for any value of 1D Fermi wave vector.

Conclusions
The purpose of this work was to show the fundamental limitation of the 3D local/semilocal exchangecorrelation energy functional approximations of DFT by considering systems with 2D characteristics. We have shown that the dimensional crossover from 3D to 2D of the exact xc energy can be significantly improved at a meta-GGA level, and we derive different exact constraints using an IBM quasi-2D. The 2D limit can be considered as a constraint on approximate functionals. For the 1D limit case we have obtained the F x ∝ 1/s 2 constraint with the IBM quasi-1D model and we have proposed a new functional that works well in this limit: the Q1D GGA functional.

Author contributions
This work was carried out with the same contribution of the two authors, in particular Lucian Constantin dealt with sections 2 and 3 in which he exhibited the quasi-2D and the quasi-1D models whose calculations and plots were made by both previously while Vittoria Urso took care of the abstract, the introduction and the conclusions.
Funding Open access funding provided by Istituto Italiano di Tecnologia within the CRUI-CARE Agreement.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The data have not been deposited and will not be deposited because we are still working on them for further developments.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.