Semilocal Exchange Functionals With Improved Performances: The Modified Enhancement Factor For Two Dimensional Quantum Systems
aa r X i v : . [ c ond - m a t . m t r l - s c i ] A p r Semilocal Exchange Functionals With Improved Performances: The Modified Enhancement Factor For TwoDimensional Quantum Systems
Subrata Jana, Abhilash Patra and Prasanjit Samal School of Physical Sciences, National Institute of Science Education and Research, Bhubaneswar 752050,Homi Bhava National Institute, INDIA (Dated: 4 March 2018)
Semilocal exchange-correlation functionals are the most accurate, realistic and widely used ones to describe thecomplex many-electron effects of two-dimensional quantum systems. Beyond local density approximation, thegeneralized gradient approximations(GGAs) are designed using reduced density gradient as main ingredient.An enhancement factor is constructed using the inhomogeneity parameter of GGAs by taking care of thelow and high density behaviors of it. Thus, the exchange energy functional proposed by making use ofthe aforementioned enhancement factor, significantly reduces the error compare to the previously proposedgradient approximations. Another enhancement factor and corresponding energy functional is also constructedusing the inhomogeneity parameter originally introduced by Becke [J. Chem. Phys. 109, 2092 (1998)].Comprehensive testing and performance of both the functionals are demonstrated with respect to the exactexchange formalism by considering two-dimensional parabolically confined quantum dots with varying particlenumber and confinement strength as a test case.
I. INTRODUCTION
In the arena of low-dimensional research, theHohenberg-Kohn-Sham variant of density-functionaltheory(DFT) is now the most widely applied formal-ism for electronic-structure calculations. The successof DFT is due to the development of several accurateapproximations for the exchange and correlation(XC)functionals . The applications of DFT are vastlyextensive because of the construction of very accuratesemilocal density functionals . Despite of its grandsuccess, the three dimensional(3D) XC-functionals inprinciple cannot be extended directly to low-dimensionalsystems due to various limitations . The developmentof new functionals in two dimensions(2D), is an ac-tive area of research with promisingly new perspectives.However, present day studies involving low-dimensionalsystems e.g. carbenoid, graphene related materials,silicon nanowire based bio-sensors and particularly semi-conductor layers and surfaces, quantum Hall devices andvarious types of quantum dots, have keenly attracted theattention of researchers and gained momentum. In thisregard, many-body effects in low dimensions need to beaddressed properly for its greater impact in solid-stateand materials research. But due to the aforementionedtime lag between the inception of the 3D and 2D XCfunctionals, the latter has not been so successful. It isonly during the last decade or so, increased attention isbeing paid in developing 2D XC functionals.The starting point of 2D XC density functional is obvi-ously the local density approximation(LDA) . The 2D-LDA for exchange combined with the 2D correlation leads to intriguing results and establishes its superiorityover quantum Monte Carlo simulations . Subsequentattempts have also been made to reduce the errors of2D-LDA . So generalized gradient approximations(2D-GGA) were the next effective attempts in thatdirection. The 2D-GGA reduces the mean percent- age error compare to 2D-LDA. The 2D-GGAs have beenconstructed by extending Becke’s proposal to the lowdimensional regime. As a matter of which, several reli-able and accurate semilocal functionals have beenconstructed. However, non of the above functional havesatisfactorily described systems both at the low as wellas high density limit.In DFT, the degree of inhomogeneity associated withthe system is included in the construction of XC function-als through the reduced density gradient( s ), which arethe main ingredient of GGA functional. In case of slowlyvarying density, reduced density gradient approachesto zero. Thus, there are two mainstream approachesfor constructing the exchange energy functionals: and . The functionals proposed by Becke, contain exchange hole potential and using it one in prin-ciple can construct the corresponding energy functionals.But, in these cases the potentials are not the functionalderivative of exchange energies. Whereas, in case of func-tionals proposed by Perdew et al are based onthe enhancement factors along with LDA for XC. Soin contrast, to Becke’s approach, the XC potentials arenothing but the functional derivative of the correspond-ing functionals. On applying spin density scaling, onecan easily construct its spin polarization version. UnlikeGGA, the meta-GGA exchange energy functionals usethe non-interacting positive definite KS kinetic energydensity ( τ ) and ’ s ’ as its ingredients. Thus ’ s ’ togetherwith ’ τ ’ forms the higher order rung of XC functionals.In stead of ’ τ ’, Becke proposed that a new inhomogeneityparameter can be used to construct the XC-functionals.The present work aims at constructing reliable and mostappropriate enhancement factors using ’ s ’ and Becke’sinhomogeneity parameter for 2D quantum systems. Sousing the above enhancement factors we have proposedtwo semilocal exchange functionals.This work is organized as follows: In the next section,we will briefly discuss exchange hole and its connection1ith exchange energy. This will be used in the follow-ing section to construct the low and high density limitof enhancement factor. Then, we will propose a form forthe enhancement factor through extrapolation betweenthe low and high density limit. To fit and test the per-formance of the newly constructed functional, it’ll be ap-plied to study few electron quantum dots. An inhomo-geneity parameter based on coordinate transformationis also proposed which can be further used to constructseries of enhancement factors and functionals. In theappendices, we’ll illustrate a scheme for constructing po-tentials for GGA and meta-GGA energy functional usedin the present work. II. EXCHANGE HOLE AND EXCHANGEENERGY
The exchange energy is considered as the electrostaticinteraction between the electron at ~r with the exchangehole at ~r + ~u surrounding the electron. So the spin-unpolarized exchange energy can be defined as E x [ ρ ] = 12 Z d r Z d u ρ ( ~r ) ρ x ( ~r, ~r + ~u ) u . (1)The exchange hole appearing in Eq.(1), is associated withthe 1 st order reduced density matrix and is given by ρ x ( ~r, ~r + ~u ) = − | Γ ( ~r, ~r + ~u ) | ρ ( ~r ) , (2)with Γ ( ~r, ~r + ~u ) = 2 P occi ψ ∗ i ( ~r ) ψ i ( ~r + ~u ), where ψ i arethe occupied KS orbitals. The exchange hole has twoimportant properties: (i) the normalization sum rule R ρ x ( ~r, ~r + ~u ) d u = − ρ x ( ~r, ~r + ~u ) ≤
0. In 2 D , the exchange energy, E x involves the cylindrical average of the exchange hole, h ρ x ( ~r, ~r + ~u ) i cyl over the direction of ~u , i.e. h ρ x ( ~r, ~r + ~u ) i cyl = Z d Ω u π ρ x ( ~r, ~r + ~u ) . (3)Using spin-scaling relation, the exchange energy func-tional can be easily generalized to any spin polarization,i.e. E x [ ρ ↑ , ρ ↓ ] = 12 E x [2 ρ ↑ ] + 12 E x [2 ρ ↓ ] . (4)The exchange energy functional also satisfies the uniformcoordinate scaling property i.e. E x [ ρ γ ] = γE x [ ρ ] , (5)where ρ γ = γ ρ is the scaling of the electronic density.Since in terms of the enhancement factor, the GGA func-tional is given by E GGAx [ ρ ] = Z d r A x ρ ( ~r ) / F x [ s ] , (6) where A x = π ) / π and reduced density gradient, s = | ~ ∇ ρ | π ) / ρ / (which is the main ingredient of GGA func-tional). Thus, the functional, F x ( s ) must reduces tounity (1) when s = 0, in order to recover the correctexchange energy for uniform density i.e. LDA. Actu-ally, there happens to be several ways of constructingenhancement factor F x : (i) it can be constructed by us-ing the small and large gradient approximations of F x and then by employing extrapolation between these twolimits, (ii) by using properties of exchange potential orexchange energy and (iii) by imposing relevant physicalconstraints. The GGA constructed by R¨as¨anan et. al. used the approach (i). Later 2D-B88 formed by apply-ing approach (ii). Here, in this case we have constructedtwo new semilocal exchange functionals by employing theapproach (i). For doing that, we’ll now elaborate on thelow and high density limits of the enhancement factor. A. Small Gradient Behavior
Lets begin with the small gradient expansion of the en-hancement factor for exchange energy in 2D. To do this,we have revisited the formalism originally proposed byBecke in 3D and the extension of it to 2D . As the Tay-lor series expansion of cylindrical averaged conventional(because no coordinate transformation is involved.) ex-change hole is h ρ x D i = − ρ ( ~r )2 − h ∇ ρ ( ~r ) − τ + 14 | ~ ∇ ρ ( ~r ) | ρ ( ~r ) i u . (7)Now, an exchange hole, based on coordinate transforma-tion can also be proposed which is given in Appendix-(A).For small inhomogeneity, one can consider 2D homoge-neous electron gas (2D-HEG) as a good reference system.Then, the cylindrical averaged uniform exchange hole isgiven by h ρ unifx D i = 2 J ( k F u ) k F u ρ ( ~r ) , (8)where u be the separation between pair of electrons and k F = (2 πρ ) is the Thomas-Fermi wavevector in 2D. Sothe cylindrical averaged exchange hole can be expressedin terms of the polynomials of u , h ρ x D i = h a ( ~r ) u + b ( ~r ) u + ........ i h ρ unifx D i . (9)Now, truncating the polynomial up to u and compar-ing it with the Taylor series expansion of the cylindricalaveraged exchange hole i.e. Eq.(7) leads to a ( ~r ) = 12 ρ h ∇ ρ + 14 | ~ ∇ ρ | ρ i . (10)2owever, by applying normalized sum rule constraint tothe cylindrical averaged exchange hole, the coefficient b turns out to be b ( ~r ) = − π I (1) I (3) ρ ( ~r ) a ( ~r ) , (11)where I ( m ) is nothing but I ( m ) = Z z dyJ ( y ) , (12)with z corresponding to the 1 st zero of the Bessel func-tion and has to be evaluated numerically. For slowlyvarying density, semi-classical approximation of kineticenergy density can be used and upon substituting it inthe enhancement factor, the same modifies to F SGLx = 1 + µ SGL s , (13)where s = |∇ ρ | k F ρ be the reduced density gradient.Whereas, µ = 3 π / κ SGL is the small gradient coeffi-cient of the enhancement factor, with κ SGL = 4˜ κ/ . κ obtained from / √ π h I (0) I (3) − I (1) I (2) I (3) i . B. Large Gradient Behavior
The large gradient behavior of exchange hole as dis-cussed by Becke’s and R¨as¨anan which give rise to h ρ x D i ≈ h | ~ ∇ ρ | ρ R i e − α ( ~r ) R . (14)This Gaussian approximation of exchange hole is pro-posed in order to produce correct short-range behaviorand finite exchange energy at large density gradient. Theparameter α , is obtained from the normalization condi-tion of exchange hole and is given by α ( ~r ) = h πG (3)2 | ~ ∇ ρ | ρ i , (15)with G ( m ) = Z ∞ dy y m e − y . (16)On using these large gradient limit results, the new en-hancement factor can be obtained as F LGLx = 1 + µ LGL s . (17)The factor µ LGL = ( π ) κ LGL A x is the large gradient coeffi-cient of the enhancement factor, with A x = π ) π . Theparameter κ LGL is obtained to be 0 . s E nh a n ce m e n t F ac t o r F x ( s ) GGAB88MODGGA
FIG. 1. Plotted are the enhancement factors F x of the MOD-GGA functional and that of 2D-B88 and 2D-GGA for com-parison. C. The Modified Enhancement Factor
It is trivial from the SGL and LGL of the enhance-ment factor that F x behaves as s and s respectively.Now we’ll combine these LGL and SGL results to findan analogous and more general expression for enhance-ment factor that interpolates between the two limits. Apossible expression is F MOD − GGAx ( s ) = 1 + µ s log( g )1 + βs log( g ) + (1 − e − cs )(18)with g = s + √ s . This form obeys the large andsmall gradient behaviors of enhancement factor, thoughthe parameters are different from its SGL and LGL valueas obtained above. The three parameters µ , β and c areobtained by using LGL of enhancement factor by con-sidering the physically relevant 2D systems like few elec-trons parabolic quantum dots. As matter of which, theparameters µ , β and c are obtained to be 0 . µ LGL ,0 .
248 and 0 . and2D-GGA . The results obtained with the new functionalare given in the Table-(I) which confirms the significantamount of reduction of error compare to the GGA func-tionals. Henceforth, we name this functional as modifiedGGA (MOD-GGA) as it is constructed by making mod-ification over the existing GGA functional.3 ABLE I. Shown are exchange energies (in atomic units) forparabolically confined few-electron quantum dots using self-consistent calculation with OCTOPUS code . The 1 st col-umn contains the number of particles. The 2 nd column con-tains different confinement strength used for fitting the pa-rameters of the new functional. For the MOD-GGA calcu-lation, the output spin polarized density and kinetic energydensity of KLI-OEP is used. The results of the MOD − GGA functional is presented at the last column. Results for 2D
EXX , LDA , GGA and B are also shown for compar-ison. The last row contains the mean percentage error, ∆.N ω − E EXXx − E LDAx − E GGAx − E B x − E MODGGAx / . r τ W σ / τ K S σ FIG. 2. Ratio of spin polarized von Weizs¨acker (WV) kineticenergy density ( τ Wσ = | ~ ∇ ρ σ | ρ σ ) to the exact kinetic energy den-sity ( τ KSσ ) for two electrons confined in a parabolic quantumdot of confinement strength ω = 1. This figure shows thatnear origin τ KSσ → τ Wσ . III. ENHANCEMENT FACTOR FROMBECKE’S INHOMOGENEITY PARAMETER
Becke showed that the coefficient of u of the Taylorseries expansion of exchange hole in Eq.(7) is a ”self-interaction” free term i.e., the interaction between theelectron and the hole surrounding it at each referencepoint is zero for one electron. For one electron, the ki-netic energy term present within the square bracket ofEq.(7) exactly cancels with the gradient term. Thus animportant inhomogeneity parameter can be given alongwith reduced density gradient, Q B , as Q B = 1 τ h τ − τ + 18 | ~ ∇ ρ | ρ + 14 ∇ ρ i (19)with, τ = 14 k F ρ . (20)This inhomogeneity parameter can be used to definethe diffuse or compact exchange hole surrounding an elec-tron. If the exact quadratic term i.e. the term contain-ing inhomogeneity parameter is larger than the homoge-neous counterpart then it represents more compact holeif not then diffuse exchange hole. The inhomogeneityparameter is zero for uniform density. Near the origin, τ KS ≈ τ W (see figure-(2)) and therefore it depends onlyon the Laplacian of density. Also in the exponential tailregion as the KS KE density equals to the VW correc-tion, due to its one electron like character, it dependsonly on the Laplacian of density. For two dimensionalquantum system , near r → ∞ . It is alsoinvariant under uniform coordinate scaling, i.e., Q B [ ρ λ ; ~r ] = Q B [ ρ ; λ~r ] . (21)In the intermediate region, it becomes positive.4 . Modeling The Enhancement Factor An enhancement factor cab be designed using the prop-erty of the Q B . For slowly varying density, Q B is small.Thus, to recover the gradient expansion of the enhance-ment factor, one may Taylor expand F x ( Q B ) as a powerseries of Q B about Q B = 0: F x ( Q B ) = F x ( Q B ) | Q B =0 + F ′ x ( Q B ) | Q B =0 Q B + ..... (22)where F ′ x (0) = dF x dQ B (cid:12)(cid:12)(cid:12) Q B =0 . (23)As for slowly varying density, the gradient expansion ofkinetic energy density is given by, τ GEA = τ + 16 ∇ ρ . (24)Therefore, Q B = 1 τ h τ − τ + 18 | ~ ∇ ρ | ρ + 14 ∇ ρ i = 1 τ h | ~ ∇ ρ | ρ + 112 ∇ ρ i = 23 h p + 2 q i , (25)where two dimensionless parameter p and q are defined asreduced density gradient and reduced Laplacian densitygradient given by p = | ~ ∇ ρ | (2 k F ρ ) ; q = ∇ ρ k F ρ . (26)Thus, F x [ p, q ] = 1 + 23 F ′ x (0)(3 p + 2 q ) (27)with the corresponding exchange energy functional givenby E x [ ρ ] = Z d rǫ LDAx F x [ p, q ] . (28)Now to eliminate the Laplacian we use integrating byparts, so that Z d rρ ( ~r ) ǫ LDAx q = 12 Z d rρ ( ~r ) ǫ LDAx p . (29)As a matter of which, the enhancement factor becomes F x [ p ] = 1 + 23 F ′ x (0)(3 p + p )= 1 + 83 pF ′ x (0) (30)For slowly varying density limit. From Eq.(13), F x [ p ] = 1 + µ SGL p . (31) Comparing Eq.(30) and Eq.(31), we have obtained F ′ x (0) = µ SGL . So the simplest conceivable enhance-ment factor is, F x [ Q B ] = 1 + αQ B p γQ B ) (32)with α = µ SGL and γ value need to be chosen so as toreduce the mean percentage error for the overall test set. B. Performance Of The Functional
To test the accuracy and efficiency of the newly con-structed functional described above, we have applied itto the few electron parabolic quantum dot. From the testset we have chosen the value of γ to be 0 . . The results areshown in Table - II, where the new functional is denotedas M GGA . TABLE II. The table caption is same as Table-I except thatcomparison of results are done w.r.t. 2D
EXX and BR .N ω − E EXXx − E BRx E MGGAx / . IV. CONCLUDING REMARKS
To summarize, we have obtained two forms of the en-hancement factor and therefore the corresponding semilo-cal exchange energy functionals using ingredient of
GGA V X , σ ( a . u . ) OEP-KLIGGAMODGGAN=2, ω=1.5 r (a.u.) -1.5-1-0.50 V X , σ ( a . u . ) OEP-KLIGGAMODGGA N=6, ω=1 V X , σ ( a . u . ) OEP-KLIGGAMODGGA N=20, ω=1.5 V X , σ ( a . u ) OEP-KLIGGAMODGGAN=56, ω=2.5 V X , σ ( a . u ) OEP-KLIGGAMODGGAN=110, ω=1 r (a.u.) -0.06-0.05-0.04-0.03-0.02-0.010 V X m e t a - GGA ( a . u . ) meta-GGA FIG. 3. First five figures represent exchange potentials v MOD − GGAxσ of Eq.(B2) compared with GGA potential and KLI-OEPfor different number of electrons and confinement strength ω in parabolic potential. Last figure represents potential of Eq.(C7)for two non-interacting electron confined in a parabolic potential. and inhomogeneity parameter defined by Becke. Thenewly constructed functionals have reduced the errors forthe overall test set i.e., parabolic quantum dot with vary-ing particle number and confinement strength comparedto existing ones i.e., 2D-GGA, 2D-B88, 2D-BR. We havealso obtained the form of exchange potential for our pro-posed functionals. The parameters of the functionals areobtained by using low and high density limits of the en-hancement factor and by fitting with the exact exchangeresults of the parabolic quantum dot. So the proposedfunctionals will enable us for precise many-electron cal-culations of larger structures such as arrays of quantumdots and quantum-Hall devices. We believe that, the con-struction takes the GGA proposed by R¨as¨anan et. al one step forward in view of the improvement in the exchangeenergy. We have also shown that using inhomogeneityparameter of Becke, semilocal density functionals withimprove performance can be constructed. In this way,one can propose hybrid density functional for 2D systemsas it has already been designed accurately in 3D. V. ACKNOWLEDGMENTS
The authors would like to acknowledge for the financialsupport from the Department of Atomic Energy, Govern-ment of India6 ppendix A: Coordinate TransformedExchange-Hole Based Enhancement factor
Since the exchange energy is related to exchange holeand exchange hole is related first order reduced densitymatrix. Thus, the different form of exchange hole den-sity plays a significant role in designing the exchangeenergy functional. It is of great interest to study theshort-range behavior of exchange hole under general co-ordinate transformation i.e., ( ~r , ~r ) → ( ~r λ , u ), where, ~r λ = λ~r + (1 − λ ) ~r . Now the exchange energy func-tional becomes E Dx = 12 Z d r λ ρ ( ~r λ ) Z d u ρ tx D ( ~r λ , u ) u (A1)where, ρ tx D is the transformed exchange hole densitydefined by ρ tx D = − | Γ D t ( ~r λ − (1 − λ ) ~u, ~r λ + λ~u ) | ρ ( ~r ) (A2)with Γ D t , be the KS single particle density matrix. Thereal parameter λ takes the value 1 for conventional ex-change hole and for on the top of the exchange hole.Now the transformed single particle KS density matrix isexpressed around u = 0 as Γ D t ( ~r, ~u ) = e ~u. [ − (1 − λ ) ~ ∇ + λ~ ∇ ] Γ D t ( ~r, ~u ) | ~u =0 = e ~u. [ − (1 − λ ) ~ ∇ + λ~ ∇ ] occ X i Ψ ∗ i ( ~r λ − (1 − λ ) ~u )Ψ i ( ~r λ + λ~u ) | ~u =0 (A3)where, ~ ∇ and ~ ∇ operate on Ψ ∗ i and Ψ i respectively.Taking the cylindrical average of Taylor series expansionof Eq.(A3) yields the correct small u behavior i.e. h ρ tx D i = − ρ ( ~r )2 − h(cid:16) λ − λ + 12 (cid:17) ∇ ρ ( ~r ) − τ + 14 (cid:16) λ − (cid:17) | ~ ∇ ρ ( ~r ) | ρ ( ~r ) i u (A4)Define a dimensionless parameter Q λB , Q λB = 1 τ h (cid:16) λ − λ + 12 (cid:17) ∇ ρ + τ − τ + 18 (2 λ − | ~ ∇ ρ | ρ i (A5)Now using Eq.(A5), Eq.(A4) can be written as, h ρ tx D i = − ρ ( ~r )2 + 12 τ (cid:16) − Q λB (cid:17) u (A6) Thus instead of Beckes inhomogeneity parameter Q B a λ -dependent inhomogeneity parameter can be used inEq.(19) that leads to a family of enhancement factors. Appendix B: Ingredients of GGA Potential
Here we have derived an explicit expression for themodified exchange potential, v MOD − GGAx . As from exactspin scaling, the spin-labeled exchange potential is givenby v D − GGAxσ = δE D − GGAx [ ρ ↑ , ρ ↓ ] δρ σ ( ~r )= δE D − GGAx [ ρ ↑ , ρ ↓ ] δρ σ ( ~r ) (cid:12)(cid:12)(cid:12) ρ ( ~r )=2 ρ σ ( ~r ) , (B1)where σ = ↑ or ↓ and ρ ( ~r ) = ρ ↑ + ρ ↓ . With MOD-GGAenhancement factor the exchange potential becomes v MOD − GGAxσ = A x / ρ σ ( ~r ) / h F x ( s σ ) i + A x / ρ σ ( ~r ) / h − s σ ( ~r ) − k F ∇ ρ |∇ ρ | + 12 k F ~ ∇ ρ σ · ~ ∇| ~ ∇ ρ σ || ~ ∇ ρ σ | i dF x ( s σ ) ds σ + A x / ρ σ ( ~r ) / h − k F ) ~ ∇ ρ · ~ ∇| ~ ∇ ρ || ~ ∇ ρ | ρ + 32 s σ i d F x ( s σ ) ds σ , (B2)where A x = π ) π and the enhancement factor and it’sderivatives w.r.t reduced density gradient are given by F x ( s σ ) = 1 + µ s σ log( g σ )1 + βs σ log( g σ ) + (1 − e − cs σ ) (B3) dF x ( s σ ) ds σ = µ log( g σ ) + s σ g σ dg σ ds σ βs σ log( g σ ) + (1 − e − cs σ ) − µ s σ log( g σ )[ β s − / σ log( g σ ) + βs / σ g σ dg σ ds σ + 2 cs σ e − cs σ ][1 + βs σ log( g σ ) + (1 − e − cs σ )] (B4)and7 F x ( s σ ) ds σ = µ g σ dg σ ds σ + s σ g σ (cid:16) dg σ ds σ (cid:17) + s σ g σ d s σ ds σ βs σ log( g σ ) + (1 − e − cs σ ) − µ h log( g σ ) + s σ g σ dg σ ds σ i [ β s − / σ log( g σ ) + βs / σ g σ dg σ ds σ + 2 cs σ e − cs σ ][1 + βs σ log( g σ ) + (1 − e − cs σ )] − µ s σ log( g σ )[ − β s − / σ log( g σ ) + βs − / σ g σ dg σ ds σ − βs / σ g σ ( dg σ ds σ ) + βg σ d g σ ds σ − ce − cs σ − c s σ e − cs σ ][1 + βs σ log( g σ ) + (1 − e − cs σ )] + 2 µ s σ log( g σ )[ β s − / σ log( g σ ) + βs / σ g σ dg σ ds σ + 2 cs σ e − cs σ ] [1 + βs σ log( g σ ) + (1 − e − cs σ )] . (B5) Appendix C: Ingredients of Two ElectronsMeta-GGA Potential
If any general density functional is given by F [ ρ ] = Z d r G [ ρ, ~ ∇ ρ, ∇ ρ, ......., ∇ m ρ ; ~r ] . (C1)Then, the functional derivative of the above functional is δF [ ρ ] δρ = ∂G [ ρ ] ∂ρ − ~ ∇ ρ · ∂G∂ ~ ∇ ρ + ~ ∇ ρ · ∂G∂ ~ ∇ ρ + ... +( − m ~ ∇ m ρ · ∂G∂ ~ ∇ m ρ . (C2)So for the exchange energy functional E x [ ρ ] = A x Z d r ρ ( ~r ) F x [ Q B ] , (C3)the corresponding exchange potential is given by v x A x = 32 ρ ( ~r ) F x [ Q B ] + ρ ( ~r ) dF x dQ B ∂Q B ∂ρ − ~ ∇ · h ρ ( ~r ) dF x dQ B ∂Q B ∂ ~ ∇ ρ i + ~ ∇ · h ρ ( ~r ) dF x dQ B ∂Q B ∂ ~ ∇ ρ i . (C4)As for two electron systems, τ = 18 | ~ ∇ ρ | ρ . (C5)Therefore, Q B = 1 + 14 τ ∇ ρ = 1 + 14 c f ∇ ρρ , (C6)where c f = (2 π ) . Now using Eq.(C6) into Eq.(C4)exchange only potential is obtained to be v x = A x h ρ ( ~r ) F x − c f ∇ ρρ dF x dQ B + 14 c f ∇ h dF x dQ B ii (C7)with dF x dQ B = α p γQ B ) − αγ Q B [1 + ( γQ B ) ] / . (C8) For the case of two non-interacting electrons confined ina parabolic quantum dot, the electron density is given by ρ ( ~r ) = 2 π exp( − r ) (C9)In cylindrical coordinate the Laplacian operator ∇ canbe written as ∇ = ∂ ∂r + 1 r ∂∂r (C10)So by making use of the above density and the Lapla-cian in Q B one can obtain v x for a non-interacting two-electron system. REFERENCES P. Hohenberg and W. Kohn, Phys. Rev. , B864(1964). W. Kohn and L.J. Sham, Phys. Rev. , A1133(1965). A.D. Becke, Int. J. Quantum Chem. , 1915 (1983). J.P. Perdew, Phys. Rev. Lett. , 1665 (1985) J.P. Perdew and Y. Wang, Phys. Rev. B , 8800(1986). A.D. Becke, Phys. Rev. A , 3098 (1988). A.D. Becke and M. R. Roussel, Phys. Rev. A , 3761(1989). A. D. Becke, J. Chem. Phys. , 1040 (1996). J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). R.M. Koehl, G.K. Odom, and G.E. Scuseria, Mol.Phys. , 835 (1996). T.V. Voorhis and G.E. Scuseria, J. Chem. Phys. ,400 (1998). F.A. Hamprecht, A.J. Cohen, D.J. Tozer, and N.C.Handy, J. Chem. Phys. , 6264 (1998). T. Tsuneda and K. Hirao, Phys. Rev. B , 15527(2000). J. Tao, J.P. Perdew, V.N. Staroverov, and G.E. Scuse-ria, Phys. Rev. Lett. , 146401 (2003). Y. Zhao and D.G. Truhlar, J. Chem. Phys. , 194101(2006). J.P. Perdew, A. Ruzsinszky, G.I. Csonka, L.A. Con-stantin, and J. Sun, Phys. Rev. Lett. , 026403(2009).8 F. Tran and P. Blaha, Phys. Rev. Lett. , 226401(2009). J. Sun, A. Ruzsinszky, and J.P. Perdew, Phys. Rev.Lett. , 036402 (2015). J. Tao and and Y. Mo, Phys. Rev. Lett. , 073001(2016). L. P. Kouwenhoven, D. G. Austing, and S. Tarucha,Rep. Prog. Phys. , 701 (2001) S. M. Reimann and M. Manninen, Rev. Mod. Phys. ,1283 (2002). Y.-H. Kim, I.-H. Lee, S. Nagaraja, J.-P. Leburton, R.Q. Hood, and R. M. Martin, Phys. Rev. B , 5202(2000). A. K. Rajagopal and J. C. Kimball, Phys. Rev. B ,2819 (1977). B. Tanatar and D. M. Ceperley, Phys. Rev. B , 5005(1989). C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. B.Bachelet, Phys. Rev. Lett. , 256601 (2002). H. Saarikoski, E. R¨as¨anen, S. Siljamki, A. Harju, M. J.Puska, and R. M. Nieminen, Phys. Rev. B , 205327(2003). S. Pittalis, E. R¨as¨anen, N. Helbig, and E. K. U. Gross,Phys. Rev. B , 235314 (2007). S. Pittalis, E. R¨as¨anen, J. G. Vilhena and M. A. L.Marques, Phys. Rev. A , 012503 (2009). S. Pittalis, E. R¨as¨anen and E. K. U. Gross, Phys. Rev.A , 032515 (2009). S. Pittalis and E. R¨as¨anen, Phys. Rev. B , 165112(2009). S. Pittalis, E. R¨as¨anen and C. R. Proetto, Phys. Rev.B , 115108 (2010). E. R¨as¨anen, S. Pittalis, Physica E , 12321235 (2010). S. Sakiroglu and E. R¨as¨anen, Phys. Rev. A , 012505(2010). S. Pittalis and E. R¨as¨anen, Phys. Rev. B , 165123(2010). E. R¨as¨anen, S. Pittalis, J. G. Vilhena, M. A. L. Mar-ques, Int. J. Quantum Chem., , 23082314 (2010). A. Putaja, E. R¨as¨anen, R. van Leeuwen, J. G. Vilhenaand M. A. L. Marques, Phys. Rev. B , 165101 (2012). J. G. Vilhena,E. R¨as¨anen, M. A. L. Marques and S. Pit-talis, J. Chem. Theory Comput. , 18371842 (2014). S. Pittalis, E. R¨as¨anen and M. A. L. Marques, Phys.Rev. B , 195322 (2008). S. Pittalis, E. R¨as¨anen, C. R. Proetto and E. K. U.Gross, Phys. Rev. B , 085316 (2009). E. R¨as¨anen, S. Pittalis and C. R. Proetto, Phys. Rev.B , 195103 (2010). A. Becke, J. Chem. Phys. , 2092 (1998). J.W. Negele and D. Vautherin, Phys. Rev. C , 1472(1972). J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A , 5453 (1992). N. Helbig, S. Kurth, S. Pittalis, E. R¨a s ¨anen, and E. K.U. Gross, Phys. Rev. B , 245106 (2008). M. A. L. Marques, A. Castro, G. F. Bertsch, and A.Rubio, Comput. Phys. Commun.151