Modified Fourth-Order Kinetic Energy Gradient ... - ACS Publications

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Modified fourth-order kinetic energy gradient expansion with Hartree potential dependent coefficients Lucian A. Constantin, Eduardo Fabiano, and Fabio Della Sala J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00705 • Publication Date (Web): 21 Aug 2017 Downloaded from http://pubs.acs.org on August 22, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Modified fourth-order kinetic energy gradient expansion with Hartree potential dependent coefficients Lucian A. Constantin,† Eduardo Fabiano,‡,† and Fabio Della Sala∗,‡,† Center for Biomolecular Nanotechnologies @UNILE, Istituto Italiano di Tecnologia, Via Barsanti, I-73010 Arnesano, Italy, and Institute for Microelectronic and Microsystems (CNR-IMM), Via Monteroni, Campus Unisalento, 73100 Lecce, Italy E-mail: [email protected]

Abstract Using the semiclassical neutral atom theory, we developed a modified fourth-order kinetic energy (KE) gradient expansion (GE4m) that keeps unchanged all the linearresponse terms of the uniform electron gas and gives a significant improvement with respect to the known semilocal functionals for both large atoms and jellium surfaces. On the other hand, GE4m is not accurate for light atoms: thus, we modified the GE4m coefficients making them dependent on a novel ingredient, the reduced Hartree potential, recently introduced in the Journal of Chemical Physics 145, 084110 (2016), in the context of exchange functionals. The resulting KE gradient expansion functional, named uGE4m, belongs to the novel class of u-meta-generalized-gradientapproximations (uMGGA) whose members depend on the conventional ingredients (i.e. ∗

To whom correspondence should be addressed Center for Biomolecular Nanotechnologies @UNILE, Istituto Italiano di Tecnologia, Via Barsanti, I73010 Arnesano, Italy ‡ Institute for Microelectronic and Microsystems (CNR-IMM), Via Monteroni, Campus Unisalento, 73100 Lecce, Italy †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 42

the reduced gradient and Laplacian of the density) as well as on the reduced Hartree potential. To test uGE4m, we defined an appropriate benchmark (including total KE and KE differences for atoms, molecules and jellium clusters) for gradient expansion functionals, i.e. including only those systems which are mainly described by a slowly-varying density regime. While most of the GGA and meta-GGA KE functionals (we tested 18 of them) are accurate for some properties and inaccurate for others, uGE4m shows a consistently good performance for all the properties considered. This represents a qualitative boost in the KE functional development and highlights the importance of the reduced Hartree potential for the construction of next-generation KE functionals.

1

Introduction

The main quantity in the Density Functional Theory (DFT) 1,2 is the ground-state electronic density n(r) which can be found solving the Euler equation 2 δTs [n] + vext (r) + δn(r)

Z

dr′

n(r′ ) δExc [n] + = µ, |r − r′ | δn(r)

(1)

where Ts [n] is the non-interacting kinetic energy (KE) functional, vext (r) is the external potential, Exc [n] is the exchange-correlation (XC) energy functional, and µ is a Lagrange R multiplier fixed from the normalization condition dr , n(r) = N , with N being the total

number of electrons.

Equation (1) represents the original orbital-free formulation of DFT (OF-DFT), 3–7 where orbitals are not used and both Ts [n] and Exc [n] must be approximated. However, because Ts [n] gives usually the dominant contribution to the ground-state energy, 2 and because of its highly non-local nature, 5,8–12 its approximation turned out to be much more difficult than for the XC part. 5,11,13,14 In fact, the development of accurate KE functionals remains one of the biggest theoretical challenges in DFT. 15

2


Page 3 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nevertheless, important progress has been achieved in the OF-DFT, and several non-local KE functionals 16–24 have been constructed, with the main goal to yield accurate binding properties in bulk solids. 16 On the other hand, the local Thomas-Fermi (TF) functional 25,26 can not bind atoms into molecules 27 and many popular generalized-gradient approximations (GGAs) KE functionals fail badly for bulk solids. 28 More recently, Laplacian-Level metaGGA functionals (LLMGGA) have been developed. 29–35 The inclusion of the Laplacian as additional ingredient in the functional form is very important: it allows to model the fourthorder gradient expansion (GE4), to distinguish nuclear regions 33 and to accurately model the KE density. 35 Nowadays, there is a large interest in semilocal KE functionals 12,32,33,35–44 and their functional derivatives (the potential and the kernel), that can be used in different contexts in addition to OF-DFT, such as: (i) Subsystem DFT, also referred as frozen density embedding (FDE) 45–48 and its timedependent extension, 7,49,50 where an electronic system is partitioned in subsystems which interact one to the other via the orbital-free non-additive potential. Several GGA and LLMGGA functionals have been proved to be accurate for FDE calculations of weaklybounded molecular complexes, 32,51 but they fail for stronger interaction between subsystems. 52–54 (ii) Quantum hydrodynamic models 55–61 where a careful treatment of the KE kernel is essential for a correct description of plasmonic systems. (iii) Information theory 62–65 and machine learning techniques. 11,14,66 This latter approach has been in fact successfully used for the construction of XC functionals, 67–70 but the KE case has been found to be much more difficult. 11,14 The development of semilocal KE functionals is therefore a very active research field. Most semilocal KE functionals are based on modifications or extensions of the second-order gradient expansion (GE2) or GE4: see Section 2 for a short review. Gradient expansions are in fact powerful theoretical tools which describe with accuracy the KE slowly-varying

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 42

density regime of an electronic system, providing an ideal starting point for the development of approximate KE models. Therefore, the study and construction of KE gradient expansion expressions is central to improve the knowledge of the properties of the KE functional and to boost the development of semilocal KE functionals. In this article, we derive a modified fourth-order gradient expansion (named GE4m), using the semiclassical atom theory, 25,26,71–75 which becomes accurate for large, neutral, non-relativistic atoms and exact for the atom with an infinite number of electrons. This is based on the KE asymptotic expansion Ts = c0 Z 7/3 + c1 Z 2 + c2 Z 5/3 + . . . ,

(2)

where Z is the number of electrons, the leading term (i.e. c0 Z 7/3 , with c0 = 0.768745) is the Thomas-Fermi one; 25,26 this term is recovered exactly by any functional which is exact for the uniform electron gas. 30 The next two terms of the expansion (i.e. c1 Z 2 and c2 Z 5/3 , with c1 = −0.5 and c2 = 0.270 ) represent the atomic inner core Scott correction 71 and the atomic quantum oscillations, 72–75 respectively. Note that Eq. (2) is very accurate even for small Z, with errors below 1% for atoms in the periodic table. 30,36,76 The connection between the semiclassical atom theory and DFT has been established in Refs. 76,77, being based on the Thomas-Fermi density scaling 76,78 and used to develop XC 36,79 and KE 30,36 functionals. The GE4m functional is very accurate for the total KE of large atoms and jellium surfaces, but it is not accurate for atoms in the periodic table (i.e. up to Z≃100). On the other hand, another modified version of GE4 (MGE4), see Ref. 30, is very good for all atoms but not accurate for jellium surfaces. Therefore it seems that functionals depending only on the conventional ingredients (the gradient and the Laplacian of the density), can not well describe different systems. We then considered an additional ingredient, the reduced Hartree-potential (ηu ) introduced in Ref. 80 (with the modified notation in place of η u ), and made the GE4m coefficients

4


Page 5 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dependent on ηu . The resulting functional (named uGE4m), as well as other 18 semilocal functionals, were tested on a newly developed benchmark suitable for the assessment of gradient-expansion-based functionals, i.e. including only atoms, molecules and jelliumclusters which are essentially characterized by the slowly-varying density regime. The results show that the uGE4m is qualitatively superior to all other functionals, underlining the significance of the GE4m expansion and the importance of having Hartree-dependent coefficients.

2

Theory

The conventional gradient expansion of the KE, derived considering small perturbations of the uniform electron gas, is known until the sixth order. 81,82 However, the sixth order KE gradient expansion diverges for any finite systems. 82 Therefore, it is rarely considered in practical applications. On the other hand, GE2 and GE4 are popular tools to study and develop KE functionals. The GE2 functional 81,83 has the enhancement factor (see Appendix for the definition) FsGE2 = 1 +

5 2 20 s + q, 27 9

(3)

where s = |∇n|/(2(3π 2 )1/3 n4/3 ) is the reduced gradient, with n being the electron density, and q = ∇2 n/(4(3π 2 )2/3 n5/3 ) is the reduced Laplacian. The last term in Eq. (3) integrates to zero for any finite system, not contributing to the KE, and has no role in determining the kinetic potential. 35 Therefore, it is usually disregarded in most applications, although it has been shown to be important for the description of the exact positive defined KE density. 35,84–86 The fourth-order gradient expansion (GE4) 81,87 has the simplified enhancement factor FsGE4 = FsGE2 + ∆ with ∆ =

5

8 2 1 2 8 4 q − s q+ s ≥0. 81 9 243


(4)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 42

This is obtained using the Green’s theorem integration on the terms that can not be written solely as products of s and q ingredients, under the assumption that the density and its gradient vanish for r → ∞, or on any appropriate surface. For the full GE4 expression, see Refs. 81,88. The GE4 expansion is not accurate to describe the asymptotic expansion of the neutralatom energy as the atomic number Z goes to infinity. Hence, a modified version of GE4 (MGE4) has been presented in Ref. 30. The MGE4 has the enhancement factor 5 2 s − 3.841∆ , 27

FsMGE4 = 1 + 1.789 ·

(5)

which gives correctly all the coefficients in Eq. (2) and thus it is very accurate for the KE of atoms in the periodic table. However, as already mentioned in the Introduction, MGE4 gives very bad jellium surface kinetic energies. 30

2.1

Modified fourth-order gradient expansion

To overcome the limitations of the MGE4, here we consider the following generalized expression:

FsGE4m [γ, ν] = FsGE2 +

8 2 q 81

− γs2 q + νs4 ,

(6)

where γ and ν are parameters to be determined. Note that in the conventional GE4 expression, the second-order term (5s2 /27) and the fourth-order term 8q 2 /81 are obtained from the linear response theory, 5 while the other two terms present in the expression of ∆ represent higher-order response terms in jellium. Thus, the FsGE4m [γ, ν] ansatz preserves the linear response terms, while the higher-order response terms are dependent on the γ and ν parameters. In order to fix γ and ν, we first consider the evaluation of jellium surface kinetic energies. In Fig. 1 we report the GE4m[γ, ν] mean absolute relative error (MARE) for jellium surface 6


Page 7 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


kinetic energies (averaged over the bulk parameters rs = 2, 4, and 6) as a function of γ and ν. Inspection of Fig. 1 shows that there exists a region of minima, which are approximately

0.48

29.41 27.29

0.36 25.17 23.05

0.24

20.93 18.82

0.12

16.70

γ

0.00

14.58 12.46

-0.12 10.34 8.22 -0.24 6.11 3.99

-0.36

1.87 -0.48 0.01 0.07 0.13 0.19 0.25 0.31 0.37 0.43 0.49

ν

Figure 1: Mean Absolute Relative Error (MARE) for jellium surfaces kinetic energies (averaged over rs =2,4,6) for the GE4m[γ, ν] functional as a function of γ and ν. The black line represents Eq. (7). described by the equation (see the black line)

γ=ν+

5 1 8 8 − =ν+ − . 81 243 9 243

(7)

Note that this equation is also obtained requiring that GE4m[γ, ν] has the same asymptotic behavior of GE4 (in the density tail asymptotic region q → s2 such that ∆ → 5/243s4 ). Therefore, Eq. (7) is not a simple fit, but it expresses an important property of gradient expansions accurate for jellium surface kinetic energies (as e.g. GE4 is). In a second step, we fix the parameter ν using the semiclassical atom theory 30,36,79 im-

7



posing the condition ν = lim ν[Z] ,

(8)

Z→∞

with

ν[Z] =

Tsexact [Z] /

Z

−

Z

dr τ

dr τ

TF

4

TF

(FsGE2 2

8 8 5 2 + q 2 − s2 q + s q) 81 81 243

(s − s q) ,

(9)

where Tsexact [Z] is exact KE for the neutral, non-relativistic atom with Z electrons and τ T F is the Thomas-Fermi KE density (see Appendix). Clearly, all quantities under the integrals in Eq. (9) depends on r and on Z. Equation (9) has been obtained using Eq. (7) in Eq. (6) and grouping terms which depends on ν. To estimate the limit in Eq. (8), we

ν

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 42

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02

Exact fit GE4 coefficient

0

0.1

0.2

0.3

0.4

0.5

-1/3

Z

Figure 2: The values ν[Z] for noble gas atoms (with 2 ≤ Z ≤ 290) and the extrapolation of ν, found from the fitting formula ν[Z] ≈ ν + aZ −1/3 + bZ −2/3 . Values for atoms in the periodic table are with Z −1/3 > 0.2.

8


Page 9 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculate ν[Z] for noble gas atoms (with 2 ≤ Z ≤ 290) and fit the values with the formula ν[Z] = ν + aZ −1/3 + bZ −2/3 , which is physically motivated by the Thomas-Fermi scaling. 77,78 Finally, we find ν = 0.1826 ,

(10)

that is approximately five times bigger than the corresponding coefficient in the conventional GE4 (i.e. 8/243 ≈ 0.033). Then, from Eq. (7) we obtain γ = 0.2608 which is about twice the GE4 value (i.e. 1/9 ≈ 0.111). The ν[Z] values for all the considered atoms, together with the fit, are reported in Fig. 2. It can be seen that the final value of ν = 0.1826 turns out to be quite larger than the corresponding value for most of the periodic table atoms. This indicates that the semilocal fourth order gradient expansion cannot be simultaneously accurate for regular and asymptotically large atoms. Indeed, GE4m significantly overestimates the KE of the periodic table atoms (Ne-Rn), as it is also shown is Tab. 1. However, GE4m is already quite accurate for a noble gas atom with 290 electrons. Table 1: (Tsapprox − Tsexact )/TsT F for different noble gas atoms Ne Ar Rn 290e-

TF GE2 GE4 -9.35 -0.81 +0.77 -7.59 -0.60 +0.62 -4.71 -0.69 -0.07 -3.14 -0.52 -0.14

GE4m +3.92 +2.73 +0.81 +0.36

In Table 2 we report the c1 and c2 coefficients for several KE functionals, calculated using the method of Ref. 30. GE4m gives the exact c1 , but the c2 value is quite different from the exact one, being rather close to the GE2 and GE4 values instead. In GE4m c1 is exact but c2 is overestimated, confirming the overestimation of the KE for atoms of the periodic table.

Finally, Tab. 3 confirms, as anticipated in Fig. 1, that GE4m is very accurate for jellium surface being definitely better than GE4. 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 42

Table 2: The coefficients c1 , and c2 of the large-Z expansion of the KE (see Eq. (2)). Exact TF GE2 GE4 GE4m

c1 -0.500 -0.661 -0.536 -0.520 -0.500

c2 0.270 0.385 0.336 0.343 0.354

Table 3: Jellium surface kinetic energies (in erg/cm2 ), for several KE approximations and several bulk parameters rs . The reference values are from self-consistent Kohn-Sham LDA calculations. Last line reports the Mean Absolute Relative Error (MARE) in %. rs 2 4 6 MARE

TF GE2 -6110.6 -5627.6 -215.2 -170.1 -26.1 -14.8 244% 120%

GE4 -5553.0 -155.4 -9.5 64%

GE4m -5492.4 -142.1 -4.4 10%

MGE4 Ref. -5532.8 -5492.7 -190.3 -139.9 -26.4 -3.4 237%

In conclusion, we found that GE4m is accurate for both large atoms and jellium surfaces, which are the most relevant systems where a gradient expansion is expected to be very accurate. Thus GE4m can be a useful tool for DFT and it represents an important step towards the development of generally applicable accurate semilocal functionals.

2.2

The uGE4m KE functional

As pointed out above, the GE4m is accurate only for large atoms, whereas it overestimates the KE of atoms in the periodic table (i.e. for Z ≤ 100). In this section, we construct a KE functional (named uGE4m) that recovers GE4m in the limit large atoms, keeps the same gradient expansion form, but it is quite accurate for all atoms. To this purpose, let us consider the enhancement factor FsuGE4m = FsGE2 +

1 8 8 2 q − [ + (ν − )f1 ]s2 q + f2 νs4 , 81 9 243

10


(11)

Page 11 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 0 ≤ f1 , f2 ≤ 1, are positive bounded functions of the density, with the following properties: (i) In the limit of large atoms (Z → ∞), then f1 → 1 and f2 → 1, such that FsuGE4m → FsGE4m ; (ii) In order to achieve good accuracy for periodic table atoms, both f1 and f2 should screen considerably the large coefficients of the GE4m. The f2 should be close to zero for light atoms, in order to diminish the large KE overestimation of the s4 -term. Instead f1 is mostly important for light atoms, where the KE of the s2 q term may be relevant for the overall accuracy of the uGE4m. To satisfy the conditions (i) and (ii) is not an easy task, because the semilocal ingredients s and q can not distinguish between different atoms, even if they can characterize several atomic regions (e.g. at the nucleus q → −∞, in the tail of the density s, q → ∞, and in the regions with slowly-varying density s, |q| ≤∼ 1). Recently, we have introduced the Hartree-dependent atomic indicator, defined 42,80 (with a notation modification, u is changed from apex to pedex) as vu = 1/(1 + ηu ) with ηu = u/[3(3/π)1/3 n1/3 ] ,

(12)

where u(r) =

Z

dr′

n(r′ ) |r − r′ |

(13)

is the Hartree potential and ηu is the reduced Hartree potential. Note that ηu and thus vu are non-local functionals of the electronic density. Note that the Hartree potential can be easily computed in plane-waves or real-space approaches (for atomic-orbitals implementation, see e.g. Refs. 89,90), and it is anyway computed in order to solve Eq. (1). The quantity vu is invariant under the uniform density scaling (as s and q), being a good atomic indicator, as shown in Fig. 2 of Ref. 80. Thus, it seems the right ingredient for building the functions f1 and f2 of Eq. (11). Note that for large atoms (under the Thomas11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 42

Fermi scaling), vu ∼ Z −2/3 such that (vu )1/2 behaves in accord with the semiclassical atom theory. Let us consider the following functions

f1 =

1 , (1 + a(vu )1/2 )1/2

f2 =

1 , 1 + (vu )1/2 + bvu

(14)

with a and b parameters which will be fixed below. By construction, both f1 and f2 recover condition (i), and moreover they correctly behave for large atoms as f1,2 → 1 + c1,2 Z −1/3 + d1,2 Z −2/3 + . . .. The function f2 shows a large screening (as requested by condition (ii)), and the function f1 varies little over a large interval of vu values. In order to find the parameters a and b, we consider: (a) the disintegration of the hydrogen atom into two neutral atoms with fractional number of electrons, which is related to the fractional density scaling. 78 We compute the kinetic disintegration energy as

M (Q) = KE(1) − KE(Q) − KE(1 − Q),

(15)

where KE(Q) denotes the KE of the H atom with fractional electron number and nuclear charge Q, and we minimize the error of the kinetic disintegration energy. These model systems are characterized by the low-density regime; consequently, they have large values of vu (vu → 1 when Q → 0). (b) the Neon isoelectronic series. We minimize the error of the KE of ions, for various values of the nuclear charge N (10 ≤ N ≤ 1000). Such positive ions have high-density regions and they have been used in Ref. 12 to investigate the KE density at the nucleus. From these calculations, we find a = 450 and b = 240, which completes the construction of the uGE4m KE functional. In Fig. 3, we show the functions f1 and f2 , for noble atoms (He-Uuo). Because of the use of the atomic indicator vu , the curves are smooth (almost monotonic), and do not intersect 12


Page 13 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


each other. The function f1 varies slowly from He (f1 ≈ 0.09 at r = 2R) to Uuo (f1 ≈ 0.2 at r = 2R). On the other hand, f2 varies strongly from He (f2 ≈ 0.06 at r = 2R) to Uuo (f2 ≈ 0.6 at r = 2R). Both f1 and f2 have the minimum value at the nucleus (r = 0), and the maximum value far in the tail (r → ∞), where vu → 0.

3

Benchmarking the gradient expansion

As the uGE4m is a gradient expansion functional, it can only be accurate in the slowlyvarying density region (i.e. for small s values). Thus, in view of its assessment, it is of utmost importance to select the proper systems to be included in the test set. To this purpose, we consider the normalized integrated s-decomposed exact KE 1 I(s) = Ts

Z

s

t(s)ds,

(16)

0

where the s-decomposed exact KE t(s) is defined by 37

t(s) =

Z

τ (r)δ(s(r) − s)dr ,

(17)

with τ being the exact KE density and s(r) the reduced gradient at point r (note that R∞ Ts = 0 t(s)ds). The value of the function I(s) at a given value of s indicates the relative contribution to the KE of regions with a reduced gradient value smaller than s. Thus,

for a system characterized by the slowly-varying density regime, I(s) will approach 1 for s / 1 ÷ 1.5. In Fig. 4a we report I(s) for several noble gas atoms, up to 290 electrons. The plot clearly shows that Helium is not a slowly-varying system, as to obtain the 95% of the total KE contributions up to s ≈ 2.5 are required. Thus, any gradient expansion should fail for Helium, unless it works by large error cancellation (e.g. as in the case of GE2, which is very accurate for Helium). On the other hand, the Neon atom can be already considered a

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

slowly-varying system, as the 95% of the total KE can be obtained for s ≤ 1.1. Moreover, as expected from the semiclassical theory, the larger the atom the smaller s values are required. In Fig. 4b we report I(s) for several jellium clusters with different rs and number of electrons. Jellium clusters with N = 2 are not slowly-varying systems, as in the case of Helium. The behavior of larger jellium clusters (e.g. for N = 40) depends on the rs values. The systems with N ≥ 40 and rs = 2 can be considered slowly-varying systems as the 95% of the total KE can be obtained for s < 1.5. On the other hand, the case rs ≥ 4 requires a proper description of the rapidly-varying region. In Fig. 4c we report I(s) for several molecules. With the exception of H2 which is a two-electron system, all the considered molecules are to a good approximation slowly-varying systems (the 95% of the total KE can be obtained for s < 1.5). In fact, due to the presence of bonding regions, molecules are more slowly-varying systems than the corresponding atoms. Using the analysis of the normalized integrated s-decomposed exact KE, we can finally construct three significative test sets, considering for each one total KE as well as KE differences. Hence, we have: • Noble gas atoms. This set includes the Ne, Ar, Rn atoms as well as the closed-shell noble atom with 290 electrons. For these atoms, we computed the total KE and the ionization KE (IKE); • Molecules. We have considered the following molecules: NH3 , CH4 , H2 O, FH, HCN, N2 , C2 H4 , H2 CO, HOOH, F2 , SiH2 , PH, PH2 , PH3 , SO2 , ClF, HCl, SH, Cl2 , OH, O2 . All these molecules have I(s) between CH4 and Cl2 , see Fig. 4c. We calculated the total as well as the atomization KE of each molecule; • Jellium spheres. We consider the total KE and ionization KE of the magic jellium clusters with electron numbers 40, 58, 92, and 106, and bulk parameter rs = 2. Jellium clusters with rs ≤ 2 also yield slowly-varying systems when N ≥ 40, but such small rs does not correspond to any physical system (i.e. rs for simple and noble metals are in 14


Page 14 of 42

Page 15 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the range 2 − 6). Thus, in summary we have six benchmarks (three total energies and three energy differences). To assess the performance of the uGE4m functional with respect to the state of the art of semilocal KE functionals we have considered in our tests 18 functionals ranging from the simple Thomas-Fermi 25,26 local functional to some LLMGGAs. These functionals are listed and briefly described in the Appendix. Note that the I(s) values are only used to construct the benchmark sets (i.e. to select the systems) whereas they are not employed anymore to assess the performance of different functionals. Finally, to be able to compare the outcomes of the different tests for the various functionals, we have normalized the error result ei of the i-th test, dividing it by the average over all the examined functionals (with the worst one removed) hei i (in a similar way as in Ref. 92). In this way, we have obtained for each test and functional the relative performance indicator (RPI) RPIi =

ei . hei i

(18)

This measures how much the functional performs better (RPI< 1) or worse (RPI> 1) than the average of all functionals for the given test. Moreover, the RPI values from different tests are homogeneous quantities that can be easily averaged to obtain an overall assessment of the functionals.

4

Computational Details

All kinetic functionals have been computed using self-consistent (SCF) KS electronic density (i.e. in a post-SCF fashion). Calculations of noble gas atoms and jellium clusters, have been performed with the numerical Engel code, 93,94 using the Perdew-Burke-Ernzerhof (PBE) 95 XC functional.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 42

Jellium surfaces have been computed with a numerical code using the local density approximation (LDA) XC functional. Molecular calculations have been performed with the TURBOMOLE program package, 96,97 using the PBE 95 XC functional and a def2-TZVPP basis set. 98,99

5

Results and Discussion

Table 4 reports for each functional, the six RPIs of the considered benchmarks. The final two columns report the average RPI and the average ranking. The latter is computed considering, for each test, an ordered list of functionals (1 for the best, 18 for the worst) and then averaging, for each functional, the position over all tests; finally, the average value is approximated to the closest integer. In the following we discuss the results for each of the six benchmarks. Table 4: Relative performance indicator (RPI, see Section 3) for different tests and functionals. If RPI> 1 (RPI

Modified Fourth-Order Kinetic Energy Gradient ... - ACS Publications

Recommend Documents