Performance of semilocal kinetic-energy functionals for Orbital-Free

Subscriber access provided by UNIV OF SOUTHERN INDIANA

Quantum Electronic Structure

Performance of semilocal kinetic-energy functionals for Orbital-Free Density Functional Theory Lucian A. Constantin, Eduardo Fabiano, and Fabio Della Sala J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00183 • Publication Date (Web): 09 Apr 2019 Downloaded from http://pubs.acs.org on April 10, 2019

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 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 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.

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 41 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

Performance of semilocal kinetic-energy functionals for Orbital-Free Density Functional Theory 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 Microelectronics and Microsystems (CNR-IMM), Via Monteroni, Campus Unisalento, 73100 Lecce, Italy. E-mail: [email protected]

Abstract We assess several generalized gradient approximations (GGAs) and Laplacian-Level meta-GGAs (LL-MGGA) kinetic energy (KE) functionals for orbital-free Density Functional Theory calculations of bulk metals and semiconductors, considering equilibrium distances, bulk moduli, total and kinetic energies, and the electron densities. We also considered the effects of the pseudopotentials, the vacancy formation energies and the bond-lengths of molecular dimers. We found that LL-MGGA KE functionals are distinctively superior to GGA functionals, showing the importance of the Laplacian of the density in the functional construction. We extended the recently developed Pauli-Gaussian Second order and Laplacian (PGSL) functional (J. Phys. Chem. Lett. ∗

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

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

9, 4385 (2018)) including high-order corrections, achieving higher transferability and accuracy than conventional non-local functionals based on the Lindhard response function.

1

Introduction

Orbital-free density functional theory (OF-DFT) 1,2 is starting to become an attractive computational method in materials science. 3–14 Its main advantage is the linear scaling computational cost, 15,16 allowing application to very large systems. Moreover, temperature effects can be naturally included in OF-DFT, resulting in an accurate free-energy functional that can be used for simulations of systems in the warm dense matter regime. 4,17–23 OF-DFT is also attracting interest in the plasmonics community, in order to describe quantum-effects beyond the classical description in metallic nanostructures. 24–27 Nevertheless, in spite of recent progress, OF-DFT still needs important improvements, being usually not as accurate as other linear scaling methods. 28–32 The main limitation in the accuracy of OF-DFT calculations traces back to the approximated kinetic energy (KE) functional Ts [n] (with n being the electron density), due to the complicated and highly non-local nature of this functional. 33–40 Accurate and broadly applicable kinetic approximations are thus still missing and the development of KE functionals is actually a rather active research field, 15,17,41–45 including applications of machine-learning techniques. 37,41,46 The most advanced KE approximations in the OF-DFT framework are the non-local (NL) KE functionals. 34,42,43,47–62 These can satisfy several important constraints such as the linear response of the non-interacting uniform electron gas (i.e. the Lindhard function in the reciprocal space 34 ) and the high- and low-momentum limits. However, they are usually quite complicated methods, depend on several (possibly system-dependent) parameters, and are, in many cases, limited to periodic codes. In fact, most of them depend on the average density in the unit cell (n0 ), and therefore, they cannot be applied to finite systems or 2


Page 2 of 41

Page 3 of 41 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


interfaces/surfaces, as n0 is well-defined in these cases. A more attracting alternative is represented by the semilocal KE functionals. 46,63–75 Semilocal functionals which are defined through a KE density which depends, at each point in space, on the electron density and its gradient, are named Generalized-GradientApproximation (GGA) functionals. In the more sophisticated Laplacian-Level meta-GGAs (LL-MGGA) functionals, the KE density depends also on the Laplacian of the density. Both GGA and LL-MGGA functionals are therefore quite simple to construct, are computationally very efficient and can be applied to finite systems. The availability of accurate semilocal KE functionals would thus be a great benefit for OF-DFT applications. For many years, semilocal approximations have been retained to be inaccurate for OF-DFT calculations. 75 Only very recent developments 44,45 have reported performance significantly close to the best state-of-the-art NL functionals, for structural and energetical ground-state properties of bulk metals and semiconductors. In this Article we consider in detail the performance of these functionals, carrying on a thorough assessment not only for solid-state systems but also for molecular systems. We propose a novel LL-MGGA functional, by considering high-order gradient and Laplacian contributions, which is even more accurate than the one in Ref. 44 and outperforming conventional non-local functionals for different properties. We also provide a rationalization for the results of the different functionals, with special attention to the importance of the Laplacian contributions and the role of the pseudopotentials. The Article is organized as follows: In Section II we present the various functionals considered; In Section III, we discuss our benchmark test-set to asses the performances of the functionals; In Section IV we discuss the results for bulk solids, vacancy formations, and molecular dimers. Finally, in Section V we outline our conclusions.

3



2

Page 4 of 41

Theory

Any semilocal KE functional can be written in the form

Tssemiloc [n] = TsW [n] + Tsθ [n] , Z θ Ts [n] = τ T F (r)Fsθ (s, q)dr ,

(1) (2)

where s = |∇n|/[2(3π 2 )1/3 n4/3 ] and q = ∇2 n/[4(3π 2 )2/3 n5/3 ] are the reduced gradient and Laplacian, respectively, τ T F = (3/10)(3π 2 )2/3 n5/3 is the Thomas-Fermi (TF) KE density, 76,77 R TsW = τ T F (r)(5/3)s2 dr is the von Weiszäcker (W) kinetic energy 78 and Tsθ [n] is the Pauli KE functional. 79,80 Exact constraints 79 indicate that

Fsθ ≥ 0 .

(3)

This condition can be easily imposed in the construction of semilocal functionals, even if, many KE functionals in literature violate it (e.g. the second-order gradient expansion and related functionals). In this work we consider several possible approximations for Fsθ : - The family of TFλW functionals 44 that are defined by the Pauli enhancement factor 5 Fsθ,T F λW (s) = 1 + (λ − 1) s2 , 3

(4)

with λ being a parameter. The simplest member of this family is obtained by setting λ = 1 and simply corresponds to the sum of the TF and W kinetic energies (i.e. the TFW functional). The TFW tends, however, to definitely overestimate the kinetic energy, despite Eq. (3) is satisfied only when λ = 1. Indeed, a much better approximation is obtained for λ = 0.6, see Ref. 44, corresponding to the TF(0.6)W functional, which gives a rather good performance for most OF-DFT calculations.

4


Page 5 of 41 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 Pauli-Gaussian (PG) KE functionals, 44 defined by 2

Fsθ,P Gµ (s) = e−µs ,

(5)

with µ ≥ 0 being a parameter. A few choices are possible for the parameter µ. In particular, we mention µ = 1, resulting in the PG1 functional that showed the best OF-DFT performance for bulk metals and semiconductors, 44 and µ = 40/27, that defines the PGS functional which performs sightly worse than PG1 but satisfies the second-order KE gradient expansion (GE2). 81 - The Luo, Karasiev, and Trickey (LKT) 45 functional, that is defined by

Fsθ,LKT (s) = 1/ cosh(1.3s) .

(6)

This has been constructed by imposing the satisfaction of the Pauli KE constraints for pseudopotential densities. This functional performs relatively well for solid-state OF-DFT applications. As it will be shown in this Article, setting µ = 0.75 in the PG functional yields an enhancement factor very close to the LKT one (the correlation coefficient between the two, in the interval s ∈ [0; 1], is 0.999); indeed, LKT and PG0.75 perform similarly for OF-DFT calculations. - A further class of KE functionals has been obtained adding to the PGS functional a simple Laplacian-level correction, i.e. defining the PGSLβ LL-MGGA functionals 44 by

Fsθ,P GSLβ (s, q) = Fsθ,P GS (s) + βq 2 ,

(7)

where β is a parameter. If the value of this parameter is fixed to β = 8/81 the corresponding PGSL8/81 functional recovers the fourth-order linear response of jellium. Another possible choice is β = 0.25 which yields the PGSL0.25 functional that gives

5



Page 6 of 41

the best overall linear response for jellium. The PGSL0.25 meta-GGA turned out to be accurate for equilibrium volumes, bulk moduli, total energies, and densities of bulk metals and semiconductors, competing with state-of-the-art non-local KE functionals. 44

2.1

The PGSLr functional

In this work we propose an additional improvement over the PGSL0.25 functional by including higher-order correction terms from the fourth-order gradient expansion (GE4). 82–85 In order to preserve the good performance of PGSL0.25 for the linear response of jellium, we consider only terms proportional to qs2 and s4 . In fact, given a semilocal functional with Taylor expansion F s ≈ a1 + a2 s 2 + a3 q + a4 s 4 + a5 s 2 q + a6 q 2 ,

(8)

the corresponding jellium linear response function is 71 " 2 4 # 9 k 9 k π2 a1 + a2 + a6 , −1/χ ≈ kF 5 2kF 5 2kF

(9)

with kF = (3π 2 n)1/3 . Therefore, the qs2 and s4 terms do not contribute to the jellium linear response. Consequently, we consider the Pauli enhancement factor

Fsθ,λσ = Fsθ,P GSL0.25 (s, q) − λqs2 + σs4 ,

(10)

where λ and σ are parameters to be fixed and we impose the condition λ2 ≤ σ such that the whole functional is still integrable for any finite system whose density decays exponentially and it has a positive Pauli enhancement factor (Fsθ,λσ ≥ 0). To fix the numerical value of the parameters we have performed a scan for λ ∈ [0 : 0.45] and σ ∈ [0 : 0.25], considering a previously presented benchmark, 44 i.e. containing 9 metals and 9 semiconductors and four properties (volume, bulk modulus, total energy at equilibrium,

6


Page 7 of 41

and density errors). In particular, we consider the Mean Absolute Relative Error (MARE) normalized to the PGSL0.25 counterpart, i.e.

rM ARE(λ, σ) =

1 X M AREpP GSL0.25λσ , 4 p M AREpP GSL0.25

(11)

where p runs over four properties and the MARE is the mean absolute relative error for the 18 systems, see Ref. 44 for details.

0.25

1.15

0.2

1.1 1.05

0.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


1 0.1 0.95

forbidden region

0.05 0

0.9 0.85

0

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

λ Figure 1: The normalized error rM ARE(λ, σ) of Eq. (11) for the PGSL0.25λσ family of functionals defined in Eq. (10). The original PGSL0.25 functional is represented by the point λ = σ = 0. The results of the scan are presented in Fig. 1. We observe that for each value of λ the results usually are worsening when σ increases. Thus, the best results are obtained close to the critical values σ = λ2 . However, for σ = λ2 the calculations are often not converging. Therefore we have considered (and reported in the figure) only the (λ, σ) points for which all calculations are converging well. The minimum error is then given by λ = 0.4, and σ = 0.2. These values are used to define a functional named PGSLr. This meta-GGA improves over the PGSL0.25 results about 10%. Other important features of this functional 7



will be discussed in the Results section.

3

Systems and Methods

To assess the performance of the various functionals, in this work we have considered an extended benchmark with the following systems: • 16 bulk metals: Al, Mg, Li, and Si in their simple-cubic (sc), face-centered-cubic (fcc), body-centered-cubic (bcc) , and hexagonal-close-packed (hcp) structures; note that Silicon is metallic in these phases. • 10 bulk semiconductors: the cubic-diamond (cd) Si, and III-V cubic zincblende (zb) semiconductors (AlP, AlAs, AlSb, GaP, GaAs, GaSb, InP, InAs, and InAs); 50 We have used both LDA and GGA bulk-derived local pseudopotentials (BLPSs), 86,87 together with the corresponding local density approximation (LDA; in the Perdew and Zunger parametrization 88 ) or the Perdew-Burke-Ernzerhof (PBE) 89 exchange-correlation functionals. The functionals have been tested among the following properties: 1. The equilibrium cell volume (V0 ): this is the main target of any DFT calculation and it can usually be reproduced quite accurately. 2. The bulk modulus (B): this describes the curvature of the potential energy surface (PES) around the minimum. It is a very important physical property, but, on the other hand, it is very difficult to be reproduced with accuracy. Note also that the bulk moduli computed using Kohn-Sham (KS) DFT deviates from experiment by 10% or more. 90 3. The total energy (E0 ) per atom at equilibrium volume: the total energy has been often considered in the assessment of functionals for OF-DFT 42,45,47,50,91 However, it 8


Page 8 of 41

Page 9 of 41 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


is not an experimental observable, and moreover it relies on huge error cancellation between approximated density and approximated kinetic energy functional. 64,92 In fact, considering that the density from OF-DFT deviates from the exact KS density, then all the energy contributions (Coulomb, kinetic, nuclear) will be different, and, usually it is only by coincidence that the error on total energy is small. Thus, no assessment of the quality of a given KE functional considering only the E0 error can be done. On the other hand, the quality of the PES can be well benchmarked considering V0 , B and E0 altogether. 4. The density error (D0 ) defined as 44 1 D0 = Ne

Z

|nKS−DF T (r) − nOF −DF T (r)|d3 r ,

(12)

where Ne is the number of electrons in the unit cell and nKS−DF T (r) and nOF −DF T (r) are the KS and OF-DFT electron densities, respectively (both are computed at the KS-DFT lattice constant). This is the most precise test to verify the quality of the OF-DFT calculations, and directly check the quality of the kinetic potential, i.e. the functional derivative of the KE functional, which, self-consistently determines the OFDFT density. However, the errors can be quite large when semilocal functionals are considered. 44 5. The total kinetic energy (K0 ) of the unit cell, calculated non-self-consistently on the KS density at the KS equilibrium structure: This error directly benchmarks the quality of the KE energy functional. 75 Note that the error on K0 does not rely on error cancellation (as for E0 ), because the density is fixed. On the other hand K0 is not an experimental observable. The present benchmark largely extends the previous one considered in Ref. 44, including more systems (Silicon), more phases (hcp), more properties (K0 ) and more pseudo potentials (both LDA and GGA, instead of LDA for semiconductors and GGA for metals). 9



Page 10 of 41

For an overall assessment of the different functionals, the different properties have to be properly weighted and summed together: we thus defined the global indicator

η=

1 X M AREp , 5 p Wp

(13)

for both metals and semiconductors, where the weights Wp are reported in Table 1. The Table 1: Weights (Wp ) for the various properties and average MARE among all semilocal functionals (but TFW) for both metals (Amet ) and semiconductors (Asem ) considering LDA pseudopotentials (PBE pseudopotentials in brackets). Aavg is the average among all the As reported in the table. p V0 B E0 D0 K0

Wp 2 10 2 5 10

Aavg Amet Asem 3.2 4.5(4.1) 2.6(1.7) 13.5 18.0(18.5) 8.8(8.7) 1.2 0.3(0.3) 2.1(2.1) 9.2 3.8(3.5) 14.1(15.3) 7.1 1.5(1.4) 15.4(10.4)

weights can be considered as an acceptable performance for the functionals: a functional with a MARE much smaller (much larger) than Wp for the property p can be considered very good (very bad). The weights Wp have been defined as follow: we first consider the average among all semilocal functionals (but TFW, which is very bad and reported in this work only for comparison) for both metals (Amet ) and semiconductors (Asem ), see Table 1 (note that semilocal functionals perform very differently for metals and semiconductors). Then we define Wp as integer numbers in the range given by Amet and Asem . Without any bias one could consider setting Wp = Aavg , where Aavg is the average among all systems and pseudopotentials reported in Tab. 1. However, as discussed above, some properties (i.e. V0 , D0 , B) are more relevant than others (E0 ,K0 ), so that we set Wp < Aavg or Wp > Aavg , respectively. For example, considering E0 , Wp can be considered in the range 0.3 − 2.1. However, E0 relies too much on error cancellation, thus its importance in the final indicator has to be reduced: instead of the averaged value of 1.2, we set Wp = 2.

10


Page 11 of 41 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


Finally, for additional benchmarking of the functionals, we considered the following properties and systems: • Vacancy formation energies of Al-fcc, Li-bcc, and Si-cd, which have been previously considered to verify the accuracy of KE functionals. 50,75,91 • Equilibrium bond lengths of 8 molecular dimers: Al2 , Si2 , Ga2 , P2 , Mg2 , In2 , Sb2 , and Li2 in the singlet state. 93 Note that the OF-DFT calculations for molecules is a challenging problem for DFT 93,94 and many semilocal and non-local KE functionals even cannot be applied to finite systems. As an additional comparison, we have considered the OF-DFT results from the non-local SM 61 and HC 47 functionals. For HC we employed the optimized universal parameters for semiconductors 47 (λ = 0.01177, and β = 0.7143). Note that these functionals do not employ system-dependent parameters, as it is for the semilocal functionals considered in this work, and thus a safer comparison can be done. Moreover, HC does not even depend on the average density n0 , thus it can be also applied to finite systems, in contrast to SM. Other interesting functionals present in literature cannot be applied in the present context, i.e. where we aim to verify their performance for a wide range of systems and without the use of system-dependent parameters. For example, the vWGTF1 and vWGTF2 functionals of Ref. 75 and the WGC functional of Ref. 49 doesn’t bind semiconductors, while the EvW-WGC functional of Ref. 50 improves over HC but contains system-dependent parameters.

3.1

Kinetic energy decomposition

When semilocal functionals are considered it is important to verify how different systems behave with respect to semilocal ingredients, i.e. s and q. Thus, we considered the s-decomposition, t[n](s), and q-decomposition, t[n](q), of the

11



TF KE (t[n]) energy so that 44,64,85 Z

+∞

Ts [n] =

ds t[n](s) Fs (s) ,

(14)

dq t[n](q) Fs (q) .

(15)

0

Z

+∞

Ts [n] = −∞

considering the exact KS density with GGA pseudopotentials. These decompositions for all the bulk systems considered are reported in Fig. 2 1

t(s)

10

10

Al, Li, Mg Semiconductors Si (bcc, fcc, hcp) Si-sc Si-cd

0

-1

10 0

0.2

0.4

s

0.8

0.6

1

1

10

t(q)

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 41

10

0

-1

10 -1

-0.8

-0.6

-0.4

-0.2

0

q

0.2

0.4

0.6

0.8

1

Figure 2: s-decomposition (upper panel) and q-decomposition (lower panel) of the TF KE for all the bulk systems. Note the log-scale for the vertical axis. Figure 2 clearly shows that all the semiconductors considered behave very similarly each other and with the relevant range of s and q quite limited, i.e. s < 1 and −0.6 < q < 0.8. On the other hand, metallic systems show very different features and structures. This means that for metals the performance of semilocal functionals can be very different for different 12


Page 13 of 41 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


metals, see also Sec. 4.1 in the following. Figure 2 also shows that metals are characterized by a reduced ranges for the semilocal variables (s < 0.65 and −0.5 < q < 0.6) with respect to semiconductors and it confirms that Si-bcc, Si-fcc, Si-hcp, and Si-sc behave similar to other metals.

3.2

Computational Details

All OF-DFT calculations have been performed with the PROFESS 3.0 (PRinceton OrbitalFree Electronic Structure Software), where we implemented all the GGA and LL-GGA presented in section 2. For each considered system we have generated reference results by considering KS 95 calculations performed using the same computational set up as that used for the corresponding OF-DFT calculations. The reference calculations have been performed using the ABINIT code. 96 In all calculations we have used a KE cutoff of 1600 eV. Equilibrium volumes and bulk moduli have been calculated using a sixth-order Birch-Murnaghan’s equation of state. 97 For the bulk modulus we have considered 11 uniformly-spaced points around the equilibrium volume within an interval of ±5%. For molecular dimers, we took KS reference results from Ref 93. OF-DFT calculations have been done using LDA BLPS and the computational setup described in Ref. 93.

4

Results and Discussion

4.1

Results for bulk systems

In Table 2 we report, for various functionals, the MARE (in %) with respect to the KS references for different properties of metals and semiconductors calculated in OF-DFT simulations using GGA pseudopotentials and GGA XC functional. In the case of metals, all semilocal functionals (but TFW reported here only for completeness and not discussed in the following) perform similarly to the best non-local KE 13



Table 2: Mean absolute relative error (MARE, in %) for the equilibrium cell volume (V0 ), bulk modulus (B), total cell energy (E0 ), density (D0 ), total cell kinetic energy (K0 ) and the global indicator η, for metals and semiconductors. All calculations used GGA BLPS and GGA XC functional. The best (worst) results from semilocal functionals (excluding TFW) are highlighted in bold style (underlined). A star indicates the best result overall. Functional

rung

TFW GGA TF(0.6)W GGA LKT GGA PG1 GGA PG0.75 GGA PGSL8/81 MGGA PGSL0.25 MGGA PGSLr MGGA SM NL HC NL TFW GGA TF(0.6)W GGA LKT GGA PG1 GGA PG0.75 GGA PGSL8/81 MGGA PGSL0.25 MGGA PGSLr MGGA SM NL HC NL

V0

B E0 D0 K0 η Metals 6.26 48.19 1.01 5.53 3.80 1.99 3.70 27.63 0.32 3.41 0.78 1.11 3.61 23.91 0.18 3.21 0.61 1.00 *3.41 16.90 0.43 3.73 0.89 0.89 3.21 0.66 1.04 3.72 25.68 *0.16 0.53 5.15 1.11 1.14 6.22 11.89 4.56 11.02 0.33 2.98 4.27 0.91 3.68 12.27 0.23 3.01 1.80 0.79 3.82 *4.68 0.19 *1.47 *0.50 *0.56 5.84 11.74 0.43 2.46 0.91 0.97 Semiconductors 8.38 49.57 5.51 24.04 25.75 3.86 1.70 *5.79 2.97 16.55 12.41 1.49 1.33 13.46 2.82 17.39 12.76 1.63 2.12 7.08 2.01 15.25 9.61 1.36 1.16 9.58 2.97 16.96 13.04 1.54 3.10 10.43 0.88 13.42 6.52 1.27 1.61 7.16 1.70 14.16 10.53 1.25 *1.11 7.14 1.48 13.41 8.16 1.10 10.63 42.87 0.76 *8.03 3.35 2.38 1.98 11.67 *0.51 8.08 *3.02 *0.87

14


Page 14 of 41

Page 15 of 41 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


functional, (i.e. SM) for V0 , E0 and K0 . Exceptions are the poor accuracy of PGSL0.25 for K0 and the poor accuracy of PGSL8/81 for E0 . For B and D0 SM is definitely superior. For bulk modulus we found that the large errors of semilocal functionals mainly come from metallic silicon (Si-fcc, Si-bcc and Si-hcp): without those systems the MAREs for all semilocal functionals are reduced by about a factor of 2. Globally, the best semilocal functional (η=0.79) for metals is PGSLr, which is superior to PGSL0.25 for V0 , E0 and K0 . The best GGA functional is PG1, with η=0.89, which has also the best overall performance for V0 . At the GGA level of theory, we also point out the excellent performance of PG0.75 and LKT for E0 (even better than SM) and K0 . When semiconductors are considered, the performance for V0 and B of most semilocal functionals is even superior to the one of the best non-local functional (in this case HC). Note that we are here discussing the performances obtained employing GGA pseudopotentials and XC functionals: in this case we found that HC is not very accurate for V0 , and B; previously, in literature the performance of HC had been verified with LDA pseudopotentials (see next section). On the other hand, all semilocal functionals are significantly worse than non-local functionals for E0 , D0 and K0 . One exception is the good accuracy of PGSL8/81 for E0 . Globally the best semilocal functional for semiconductors is again PGSLr (η=1.10), which is superior to PGSL0.25 for all properties. At the GGA level, the best functional is PG1, with η=1.36, whereas LKT shows the worst performance. Concerning the large errors of D0 for all semilocal functionals (MARE in the range 13%16%) we point out that also the best non-local KE functional (HC) has quite large errors (MARE about 8%), so that no significant difference occurs when density plot are analyzed. In fact, as shown in Fig. 3), where we report the densities for two different semiconductors, while HC is better than semilocal functionals near the nuclei, PGSLr and PGSL0.25 are better near the bonding regions. An important feature to investigate concerns the similarities in the trends of different functionals for various properties and classes of systems. In particular, we focus on the

15



Electron density (a.u.)

0.12 0.1 0.08 KSDFT

0.06

PGSL0.25 PGSLr HC

0.04 0.02

Ga

0 0

As 0.2

Ga 0.4

0.6

0.8

1

GaAs [111] direction 0.09

Electron density (a.u.)

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 41

0.08 0.07

KSDFT PGSL0.25 PGSLr HC

0.06 0.05 0.04 0.03 0.02 0.01

Si

Si

Si

0 0

0.2

0.4

0.6

0.8

1

Si [111] direction Figure 3: Reference KS-DFT and OF-DFT self-consistent electron densities for GaAs (top) and Si-cd (bottom). 16


Page 17 of 41 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


PGSLr functional and how its errors correlate to those obtained using other KE functionals. This information is summarized in Fig. 4 where we report, for the different properties and both the metals and the semiconductors, the values of the (Pearson) correlation coefficient between PGSLr and the other functionals. In the case of semiconductors, one can readily

Volume

Bulk Modulus Metals

Semic.

Energy

Density

Kinetic Energy

Average

Figure 4: Pearson correlation coefficient between PGSLr and other functionals for various properties in the cases of bulk simple metals and semiconductors. observe that all functionals display a high degree of linear correlation for the E0 , K0 and

17



D0 . This means that they all show similar trends for these properties. This is not surprising considering the decomposition reported in Fig. 2. However, because the MAEs and MAREs are rather different, our finding suggests that the various functionals incorporate a systematic error for these properties. A slightly different situation is found for the structural properties (V0 and B) of semiconductors. In these cases, in fact, most semilocal functionals are still highly linearly correlated, but the non-local KE and TFW display different behaviors. For the latter this is due to the fact that TFW, despite underestimating significantly all the volumes and overestimating all the bulk moduli, has a more pronounced tendency to do so for the Sb compounds than for the phosphides. Similarly, for the HC functional, and partially for SM, the results of AlP, GaP, and InP display significantly larger errors than for other materials. This is the origin of a different behavior with respect to the semilocal KE functionals, where the errors for the various materials are more similar; this is also the principal cause of the quite large MARE displayed by these functionals for the structural properties of semiconductors. For metals the situation is quite different. In fact, PGSLr correlates well only with PG1, PGSL8/81, and PGSL0.25. While the latter similarity is quite understandable, because PGSLr is built on top of PGSL0.25, the very good linear correlation with PG1 and PGSL8/81 is less obvious, considering the different performance. This can be explained observing a systematic difference between PGSLr and PGSL8/81 (and PG1). Concerning the GGA and non-local functionals, we note that the linear correlation with PGSLr is quite poor (negative in many cases). In fact, as shown in Fig. 2, different metals behaves very differently, thus different functionals can show alternating good and less good results depending on the particular metal.

4.2

Comparison with LDA pseudopotentials

Similar overall results as those discussed above are obtained if the OF-DFT calculations are performed using LDA pseudopotentials, instead of GGA one. Results are collected in Table 18


Page 18 of 41

Page 19 of 41 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


Table 3: Same as Table 2, but using in OF-DFT calculations LDA BLPS and the LDA XC functional. Functional

rung

TFW GGA TF(0.6)W GGA LKT GGA PG1 GGA PG0.75 GGA PGSL8/81 MGGA PGSL0.25 MGGA PGSLr MGGA SM NL HC NL TFW GGA TF(0.6)W GGA LKT GGA PG1 GGA PG0.75 GGA PGSL8/81 MGGA PGSL0.25 MGGA PGSLr MGGA SM HC

V0

B E0 D0 K0 η Metals 7.02 45.67 1.05 5.58 3.56 2.02 4.09 26.48 0.33 3.71 1.02 1.14 4.06 22.31 0.18 3.46 0.90 1.03 *3.86 15.95 0.45 3.96 1.19 0.93 3.47 0.87 1.07 4.14 24.21 *0.17 6.43 13.16 0.53 5.25 0.92 1.19 4.84 11.69 0.36 3.21 4.00 0.96 3.90 11.89 0.25 3.23 1.71 0.82 4.06 *4.95 0.20 *1.95 *0.73 *0.62 5.95 13.18 0.46 2.81 1.15 1.04 Semiconductors 5.73 28.93 5.28 22.96 31.41 3.23 2.62 7.69 2.94 15.59 17.49 1.68 2.02 13.41 2.74 15.87 17.86 1.74 3.03 8.83 1.99 14.19 14.57 1.54 1.99 4.79 2.90 15.92 18.15 1.58 3.49 13.79 0.85 12.23 11.33 1.43 2.55 7.34 1.60 12.85 15.50 1.39 2.30 5.73 1.40 12.23 13.04 1.23 7.49 38.43 0.78 7.96 8.00 2.07 *1.55 *4.22 *0.44 *7.79 *7.66 *0.75

19



3 where the MAREs of these latter calculations are reported. The most significant difference (i.e. above 20% of deviation) are related i) to K0 for all systems and ii) to V0 for semiconductors: i) For most functionals the MARE for K0 is much larger with LDA pseudopotentials. Only exceptions are the functionals PGSL8/81, PGSL0.25, PGSLr for metals. ii) For all semilocal functionals, V0 is much better with GGA pseudopotentials, while the opposite is true for HC. One possbile explanantion can be related to the fact that parameters for HC have been optimized considering semiconductors with LDA pseudopotentials. For HC also the bulk moduli are much improved with LDA pseudopotentials. In general, for metals, results are (slightly) better using GGA pseudopotential than LDA: and this is true for all functionals considered. This is clearly evinced in Fig. 5 where we summarize the η values for metals and semiconductors for both pseudopotentials. On the other hand, for semiconductors, results are better using GGA pseudopotentials instead of LDA ones, only for semilocal functionals, whereas the opposite is true for HC and SM (see above). In general the errors for semiconductors change more significantly than for metals. We need to remark, however, that an important difference exists between performing LDA or GGA OF-DFT calculations (with the setup considered in this work) for semiconductors. In fact, despite both approaches yield a similar approximation to the reference KS results, only the GGA (i.e. PBE) KS results are quite close to the exact experimental values, whereas the LDA KS results yield strongly compressed lattice constants, in particular for semiconductors with Gallium and Indium, as shown in Fig. 6. This means that only the GGA OF-DFT calculations are suitable for practical applications, while the LDA ones must be considered with caution for practical purposes. On the other hand, with the aim of assessing the KE functionals, the use of the LDA or GGA OF-DFT scheme is basically equivalent, and a KE functional with broad applicability should work with both of them. Considering that LDA lattice constants are smaller than the GGA ones, it is worth to analyze if the different behavior of KE functionals discussed in Fig. 5 can be traced back to this. For a better understanding, a comparison LDA and GGA pseudopotentials for each

20


Page 20 of 41

Page 21 of 41

2.5 GGA LDA

SM

2.0

η Semiconductors

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


LKT TF(0.6)W

1.5

PG0.75

PG1

PGSL8/81 PGSLr

PGSL0.25

1.0 HC

1.0

0.5

1.5

η Metals Figure 5: η values for metals vs semiconductors, for all functionals considered. The η values are reported for both pseudopotentials: the arrows connect results obtained with GGA pseudopotentials to results obtained with LDA pseudopotentials. 21



6.8 BLPS-LDA BLPS-GGA VASP-LDA VASP-PBE Exp.

6.6 Lattice constant [Å]

6.4 6.2 6 5.8 5.6

InP

InAs

InSb

GaP

GaAs

GaSb

AlP

AlAs

5.2

AlSb

5.4

Si

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 22 of 41

Figure 6: Lattice constants for semiconductors from KS calculations using BLPS-LDA and BLP-GGA, VASP calculations 90 and experimental results.

22


Page 23 of 41 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


different property is required (the η values reported in Fig. 5, instead, correspond to the averaged value among five different properties): this is reported in Fig. S1 in the Supporting Information. Figure S1 shows that K0 is always better for GGA pseudopotentials while D0 is always better for LDA ones, and this is true for all functionals. As K0 and D0 directly reflect the accuracy of the total kinetic energy and the kinetic potential, respectively, we conclude that when the lattice constant is reduced (i.e. going from GGA to LDA pseudopotentials) the kinetic potential is better and total kinetic energy is worse, and this is true for semilocal as well as non-local KE functionals, which thus behave in a similar man ner. The differences in the η values reported in Fig. 5 are instead dominated by the errors in V0 (and B0 ), that depend on a delicate balance between energy and potential errors.

4.3

Results for vacancies

As a further test for the various functionals, we have considered vacancy formation energies for Al-fcc, Li-bcc, and Si-cd (see Table 4). In this case, in addition to SM and HC we show, for comparison, also the results obtained using the KGAP non-local KE functional. 43 This functional, which depends on the exact (fundamental) band gap energy (Eg ) of the material, is in fact particularly accurate for this test. In the present calculations we have used Eg = 1.17 eV for Si-cd, while we have assumed Eg = 0 for the two metals (Al-fcc and Li-bcc). Inspection of the table shows that for Li-bcc all the semilocal functionals perform quite well. This traces back to the fact that the Li-bcc density varies slowly, and the reduced gradient is always smaller than 0.5, 44 then all functionals recovering the TF limit work fairly well. However, the PGSL8/81, PGSL0.25, and PGSLr functionals, that also recover the second-order terms, yield slightly better results than others, proving the importance of this constraint. A similar situation is found for Al-fcc, but in this case the density is less slowly-varying. Thus, all GGA functionals generally overestimate the vacancy formation energy and only 23



Table 4: Errors (OF-DFT - KS-DFT, in eV) for vacancy formation energies for Al-fcc, Li-bcc, and Si-cd. The last two columns show the mean absolute error (MAE, in eV) and the MARE (in %). The best (worst) results from semilocal functionals (excluding TFW) are highlighted in bold style (underlined). A star indicates the best result overall. The last row of every panel reports the KS reference values. Functional

rung

Al-fcc Li-bcc Si-cd MAE MARE GGA pseudopotential TFW GGA 2.56 0.17 0.41 1.05 139.9 0.12 -0.96 0.88 92.3 TF(0.6)W GGA 1.55 LKT GGA 1.32 0.09 -0.88 0.76 78.8 PG1 GGA 0.98 0.07 -1.77 0.94 70.1 PG0.75 GGA 1.43 0.10 -1.02 0.85 86.4 PGSL8/81 MGGA 0.32 0.00 -1.43 0.59 30.7 PGSL0.25 MGGA 0.46 0.00 -1.36 0.60 36.8 PGSLr MGGA 0.39 0.02 -1.20 0.54 32.4 SM NL 0.02 0.04 -2.12 0.73 23.4 KGAP NL 0.02 0.04 *-0.27 *0.11 *4.9 HC NL -1.30 -0.03 -1.64 0.99 83.3 KS-DFT 0.655 0.990 3.345 LDA pseudopotential TFW GGA 2.51 0.15 *0.49 1.05 117.3 0.11 -0.86 0.83 76.8 TF(0.6)W GGA 1.51 LKT GGA 1.27 0.08 -0.74 0.70 64.5 57.9 PG1 GGA 0.94 0.07 -1.48 0.83 PG0.75 GGA 1.39 0.10 -0.81 0.77 70.8 PGSL8/81 MGGA 0.25 *0.00 -1.36 0.54 24.9 PGSL0.25 MGGA 0.39 *0.00 -1.29 0.56 30.2 PGSLr MGGA 0.33 0.01 -1.12 0.49 26.2 SM NL *0.06 0.05 -1.95 0.69 24.7 KGAP NL *0.06 0.05 *-0.18 *0.10 *6.3 HC NL -1.24 *0.00 -1.56 0.94 69.3 KS-DFT 0.782 1.018 3.200 -

24


Page 24 of 41

Page 25 of 41 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


Laplacian-level meta-GGAs perform quite accurately. In particular, the PGSL8/81 functional, that recovers the fourth-order linear response of the jellium, is the best in this case. The most difficult case is, finally, Si-cd, where the slowly-varying density regime is not very relevant. Indeed for this case also the non-local HC and SM functionals perform rather poorly and only the KGAP functional yields an accurate result. The semilocal functionals all tend to underestimate quite consistently the formation energy, with the exception of TFW that works very well, probably due to a large error compensation. Overall, the best performance among semilocal functionals is obtained by the PGSL8/81 functional, closely followed by PGSLr and PGSL0.25. The performances of these Laplacianlevel functionals are comparable to that of SM and 2-3 times better than HC.

4.4

Results for dimers

To conclude our assessment we consider in Table 5 the equilibrium bond lengths of several dimers. Note that this is a very difficult test for KE functionals because for this problem Table 5: Errors for the equilibrium bond lengths (OF-DFT - KSDFT, in ˚ A) of several dimers (LDA level). The last two columns show the mean absolute error (MAE, in ˚ A) and the MARE (in %). The best (worst) results from semilocal functionals (excluding TFW) are highlighted in bold style (underlined). A star indicates the best result overall. The last row reports the KS reference values. Note that HC results are from Ref. 93 and PGSLr values 2 are with the qs2 term damped by e−0.05s (see text). Functional Al2 Si2 Ga2 In2 P2 Sb2 Li2 Mg2 MAE MARE TFW 0.23 0.15 0.24 0.29 0.34 0.38 0.60 -0.32 *0.32 12.71 TF(0.6)W 0.17 0.13 0.19 0.22 0.35 0.42 0.39 -0.45 0.29 11.49 LKT 0.15 0.11 0.15 0.20 0.32 0.37 0.40 -0.45 0.27 10.61 PG1 0.14 0.10 0.14 0.18 0.33 0.39 0.41 -0.48 0.27 10.52 PG0.75 0.16 0.11 0.19 0.21 0.33 0.38 0.43 *-0.44 0.28 11.05 PGSL8/81 0.08 0.10 0.09 0.06 0.30 0.39 0.40 -0.56 0.25 9.51 PGSL0.25 0.11 0.04 0.06 0.10 0.29 0.35 0.42 -0.50 0.23 8.94 PGSLr 0.05 *0.03 0.07 0.09 0.28 0.29 0.43 -0.54 0.22 8.34 HC *0.02 -0.04 *0.01 *0.02 *0.20 *0.28 *0.34 -0.66 *0.20 *7.17 KS-DFT 2.473 2.284 2.323 2.644 1.942 2.431 2.781 3.405 -

the density is not fully in the slowly-varying regime. This fact may cause problems to 25



functionals based on gradient and/or polynomial expansions, that otherwise work well for bulk solids. In particular, we found that for the PGSLr functional the SCF convergence of the calculations is quite problematic because of the behavior of the −0.4qs2 term. In order to solve this numerical issue, we added a dumping factor to this term such that it becomes 2

−0.4qs2 e−0.05s . Note that this dumping procedure keeps almost unchanged the PGSLr functional in solids, where the reduced gradients are small, but helps to dump oscillations in the vacuum improving the convergence for dimer calculations. Finally we remark that for finite systems, like the dimers considered here, non-local functionals with a densityindependent kernel, e.g. SM and KGAP, cannot be applied since they depend on the average density on the unit cell, a parameter that is not well defined for finite systems. Inspection of the Table shows that all the considered GGA functionals perform rather similarly, yielding a MARE of about 11%, with PG1 and LKT giving slightly better results than the other functionals. Note that this result is almost twice as better as that (not reported) obtained using simple gradient expansion functionals, such as GE2 (MARE=21%) and the modified GE2 (MGE2, MARE=18%), which has been derived from the semiclassical atom theory. 63,98–104 A small but clear improvement is obtained moving to the Laplacian-meta-GGA level of theory. In this case in fact the functionals of the PGSL family give MAREs in the range 8-9%, with PGSLr being the best with MARE=8.3%. This result is not as good, but sufficiently close, to the one of the non-local HC functional, that displays a MARE of 7%. We recall anyway that the HC results have been obtained using optimal λ and β parameters for every dimer, instead of using an universal set of parameters. 93

4.5

Global Performances

In order to estimate the global performances of the KE functionals for all properties and systems, we first define an indicator for vacancy formation, η vac = MAREvac /W vac , and for dimer bond-lengths, η dim = MAREdim /W dim , where W vac = 50% and W dim = 10% 26


Page 26 of 41

Page 27 of 41

correspond to the averaged MARE among all semilocal functionals (but TFW). Then we define the averaged indicator as

η avg = η met + η sem + η vac + η dim

(16)

Note that η met ,η sem ,η vac in Eq. (16) have been also averaged among LDA and GGA pseudopotentials. All the η values are reported in Fig. 7.

2.5 2.35 Metals Semiconductors Vacancies Dimers Average

2 1.5

1.39

1.32

1.30

η

1.17 1.01

1

0.96

1.01 0.85

SM

HC

Lr PGS

L0.2 5 PGS

L8/8 1 PGS

.75 PG0

PG1

LKT

TFW

0

0.6

0.5

TFW

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


Figure 7: Performance indicator η for all systems and functionals. At the GGA level, the most accurate functional is the recently developed PG1 η avg = 1.17 that performs a bit better than the LKT functional, η avg = 1.30, which is found to be very close in accuracy to PG0.75. The best MGGA functional is the PGSLr with η avg = 0.85, which outperforms PGSL0.25 (η avg = 0.96). These two functionals are superior to HC, η avg = 1.01, that for vacancies is as inaccurate as GGA functionals. Note that for SM, the 27



average indicator in Eq. (16) cannot be defined, as SM cannot be applied to finite systems.

5

Conclusions

We have assessed the new generation of semilocal KE functionals which are accurate for OF-DFT calculations of bulk solids (metals and semiconductors). We considered several important properties: equilibrium volume, bulk modulus, total energy, kinetic energy, electron density, and vacancy formation energy, considering also the dependence of the results on to the pseudopotentials. Our statistical analysis of the influence of the pseudopotentials shows that all semilocal functionals behave quite differently from SM and HC. Going from the more accurate GGA pseudopotentials to less accurate LDA ones, the performance for semiconductors decreases for semilocal functionals, while it decreases for non-local ones. In addition, we have tested the semilocal KE functionals for the challenging problem of finite systems, considering the equilibrium bond-length of eight diatomic molecules in the singlet state. At GGA level, the best KE functional is PG1 of Ref. 44. The PG1 functional is better than the LKT functional of Ref. 45 which behaves very similar to PG0.75. All these GGA functionals have a positive Pauli KE enhancement factor (like TFW), but are not correct for the linear response of jellium and they are failing for vacancy formation energies (MARE over 60%). We then considered the more advanced class of Laplacian-level meta-GGA KE functionals, which includes the Laplacian of the density as an additional ingredient. The PGSLβ class of functionals achieves a much higher accuracy than GGAs for a large palette of systems and properties, including vacancy formation energies. The best member of the PGSLβ family of functionals is PGSL0.25 of Ref. 44 It satisfies important physical conditions: it is very accurate for the linear response of jellium 44 and it has a positive Pauli KE enhancement factor. The global performance of the PGSL0.25 functionals is almost equivalent to

28


Page 28 of 41

Page 29 of 41 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 non-local HC functional, which is quite inaccurate for vacancy formation. One limitation of the PGSL0.25 functional concerns the total KE which is overestimated. By adding to PGSL0.25 the contribution of the fourth-order terms qs2 and s4 , we have constructed the PGSLr functional, which preserves the accurate PGSL0.25 linear response in jellium, corrects the KE overestimation and gives an overall and systematic improvement of about 10% for all properties and systems. Indeed, the PGSLr provides the best overall results with η avg = 0.85, outperforming the non-local HC functional. These results suggest that PGSLr can be used in OF-DFT applications of large-scale systems.

Supporting Information Available Graphical representation of the difference between results using LDA and GGA pseudopotentials, for each property.

This material is available free of charge via the Internet at

http://pubs.acs.org/.

References (1) Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. 1964, 136, B864. (2) Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. Proc. Nat. Acad. Sc. 1979, 76, 6062–6065. (3) Hung, L.; Carter, E. A. Accurate simulations of metals at the mesoscale: Explicit treatment of 1 million atoms with quantum mechanics. Chem. Phys. Lett. 2009, 475, 163–170. (4) Lambert, F.; Clérouin, J.; Mazevet, S. Structural and dynamical properties of hot dense matter by a Thomas-Fermi-Dirac molecular dynamics. Europhys. Lett. 2006, 75, 681.

29



(5) Chen, M.; Hung, L.; Huang, C.; Xia, J.; Carter, E. A. The melting point of lithium: an orbital-free first-principles molecular dynamics study. Mol. Phys. 2013, 111, 3448– 3456. (6) Gavini, V.; Bhattacharya, K.; Ortiz, M. Quasi-continuum orbital-free densityfunctional theory: A route to multi-million atom non-periodic DFT calculation. J. Mech. Phys. Sol. 2007, 55, 697–718. (7) Radhakrishnan, B.; Gavini, V. Effect of cell size on the energetics of vacancies in aluminum studied via orbital-free density functional theory. Phys. Rev. B 2010, 82, 094117. (8) Gavini, V. Role of macroscopic deformations in energetics of vacancies in aluminum. Phys. Rev. Lett. 2008, 101, 205503. (9) Radhakrishnan, B.; Gavini, V. Orbital-free density functional theory study of the energetics of vacancy clustering and prismatic dislocation loop nucleation in aluminium. Phil. Mag. 2016, 96, 2468–2487. (10) Caspersen, K. J.; Carter, E. A. Finding transition states for crystalline solid–solid phase transformations. Proc. Nat. Acad. Sc. 2005, 102, 6738–6743. (11) Xiang, H.; Zhang, M.; Zhang, X.; Lu, G. Understanding Quantum Plasmonics from Time-Dependent Orbital-Free Density Functional Theory. J. Phys. Chem. C 2016, 120, 14330–14336. (12) del Rio, B. G.; Chen, M.; González, L. E.; Carter, E. A. Orbital-free density functional theory simulation of collective dynamics coupling in liquid Sn. J. Chem. Phys. 2018, 149, 094504. (13) Ho, G.; Ong, M. T.; Caspersen, K. J.; Carter, E. A. Energetics and kinetics of vacancy

30


Page 30 of 41

Page 31 of 41 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


diffusion and aggregation in shocked aluminium via orbital-free density functional theory. Phys. Chem. Chem. Phys. 2007, 9, 4951–4966. (14) Shin, I.; Carter, E. A. First-principles simulations of plasticity in body-centered-cubic magnesium–lithium alloys. Acta Mater. 2014, 64, 198–207. (15) Witt, W. C.; Beatriz, G.; Dieterich, J. M.; Carter, E. A. Orbital-free density functional theory for materials research. J. Mater. Res. 2018, 33, 777–795. (16) Carter, E. A. Challenges in modeling materials properties without experimental input. Science 2008, 321, 800–803. (17) Karasiev, V. V.; Chakraborty, D.; Trickey, S. B. Many-Electron Approaches in Physics, Chemistry and Mathematics; Springer, 2014; pp 113–134. (18) Karasiev, V. V.; Sjostrom, T.; Trickey, S. B. Finite-temperature orbital-free DFT molecular dynamics: Coupling Profess and Quantum Espresso. Comput. Phys. Comm. 2014, 185, 3240–3249. (19) Karasiev, V. V.; Chakraborty, D.; Shukruto, O. A.; Trickey, S. B. Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations. Phys. Rev. B 2013, 88, 161108. (20) Karasiev, V. V.; Sjostrom, T.; Dufty, J.; Trickey, S. Accurate homogeneous electron gas exchange-correlation free energy for local spin-density calculations. Phys. Rev. Lett. 2014, 112, 076403. (21) Karasiev, V. V.; Dufty, J. W.; Trickey, S. Nonempirical Semilocal Free-Energy Density Functional for Matter under Extreme Conditions. Phys. Rev. Lett. 2018, 120, 076401. (22) White, T.; Richardson, S.; Crowley, B.; Pattison, L.; Harris, J.; Gregori, G. Orbitalfree density-functional theory simulations of the dynamic structure factor of warm dense aluminum. Phys. Rev. Lett. 2013, 111, 175002. 31



(23) Karasiev, V. V.; Calder´ın, L.; Trickey, S. Importance of finite-temperature exchange correlation for warm dense matter calculations. Phys. Rev. E 2016, 93, 063207. (24) Xiang, H.; Zhang, X.; Neuhauser, D.; Lu, G. Size-Dependent Plasmonic Resonances from Large-Scale Quantum Simulations. J. Phys. Chem. Lett. 2014, 5, 1163–1169. (25) Toscano, G.; Straubel, J.; Kwiatkowski, A.; Rockstuhl, C.; Evers, F.; Xu, H.; Asger Mortensen, N.; Wubs, M. Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics. Nat. Comm. 2015, 6, 7132, Article. (26) Cirac`ı, C.; Della Sala, F. Quantum hydrodynamic theory for plasmonics: Impact of the electron density tail. Phys. Rev. B 2016, 93, 205405. (27) Domps, A.; Reinhard, P.-G.; Suraud, E. Time-dependent Thomas-Fermi approach for electron dynamics in metal clusters. Phys. Rev. Lett. 1998, 80, 5520. (28) Goedecker, S. Linear scaling electronic structure methods. Rev. Mod. Phys. 1999, 71, 1085. (29) Goedecker, S.; Scuseria, G. E. Linear scaling electronic structure methods in chemistry and physics. Comp. Sci. Eng. 2003, 5, 14–21. (30) Vogl, P.; Hjalmarson, H. P.; Dow, J. D. A semi-empirical tight-binding theory of the electronic structure of semiconductors. J. Phys. Chem. Sol. 1983, 44, 365–378. (31) Murrell, J. N. Semi-empirical electronic structure methods. J. Mol. Str.: THEOCHEM 1998, 424, 93–99. (32) Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; van Lenthe, E.; Witek, H. A.; Heine, T. DFTB parameters for the periodic table: Part 1, Electronic structure. J. Chem. Theory Comput. 2013, 9, 4006–4017. (33) Garc´ıa-González, P.; Alvarellos, J. E.; Chacón, E. Kinetic-energy density functional: Atoms and shell structure. Phys. Rev. A 1996, 54, 1897. 32


Page 32 of 41

Page 33 of 41 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


(34) Wang, Y. A.; Carter, E. A. Theoretical methods in condensed phase chemistry; Springer, 2002; pp 117–184. (35) Howard, I. A.; March, N. H.; Van Doren, V. E. r-and p-space electron densities and related kinetic and exchange energies in terms of s states alone for the leading term in the 1/Z expansion for nonrelativistic closed-shell atomic ions. Phys. Rev. A 2001, 63, 062501. (36) March, N.; Santamaria, R. Non-local relation between kinetic and exchange energy densities in Hartree–Fock theory. Int. J. Quantum Chem. 1991, 39, 585–592. (37) Snyder, J. C.; Rupp, M.; Hansen, K.; M¨ uller, K.-R.; Burke, K. Finding density functionals with machine learning. Phys. Rev. Lett. 2012, 108, 253002. (38) Della Sala, F.; Fabiano, E.; Constantin, L. A. Kohn-Sham kinetic energy density in the nuclear and asymptotic regions: Deviations from the von Weizsäcker behavior and applications to density functionals. Phys. Rev. B 2015, 91, 035126. (39) Constantin, L. A.; Fabiano, E.; Della Sala, F. Kinetic and Exchange Energy Densities near the Nucleus. Computation 2016, 4, 19. (40) Constantin, L. A.; Fabiano, E.; Pitarke, J.; Della Sala, F. Semilocal density functional theory with correct surface asymptotics. Phys. Rev. B 2016, 93, 115127. (41) Yao, K.; Parkhill, J. Kinetic Energy of Hydrocarbons as a Function of Electron Density and Convolutional Neural Networks. J. Chem. Theory Comput. 2016, 12, 1139–1147. (42) Mi, W.; Genova, A.; Pavanello, M. Nonlocal kinetic energy functionals by functional integration. J. Chem. Phys. 2018, 148, 184107. (43) Constantin, L. A.; Fabiano, E.; Della Sala, F. Nonlocal kinetic energy functional from the jellium-with-gap model: Applications to orbital-free density functional theory. Phys. Rev. B 2018, 97, 205137. 33



(44) Constantin, L. A.; Fabiano, E.; Della Sala, F. Semilocal Pauli–Gaussian Kinetic Functionals for Orbital-Free Density Functional Theory Calculations of Solids. J. Phys. Chem. Lett. 2018, 9, 4385–4390. (45) Luo, K.; Karasiev, V. V.; Trickey, S. A simple generalized gradient approximation for the noninteracting kinetic energy density functional. Phys. Rev. B 2018, 98, 041111. (46) Seino, J.; Kageyama, R.; Fujinami, M.; Ikabata, Y.; Nakai, H. Semi-local machinelearned kinetic energy density functional with third-order gradients of electron density. J. Chem. Phys. 2018, 148, 241705. (47) Huang, C.; Carter, E. A. Nonlocal orbital-free kinetic energy density functional for semiconductors. Phys. Rev. B 2010, 81, 045206. (48) Wang, Y. A.; Govind, N.; Carter, E. A. Orbital-free kinetic-energy functionals for the nearly free electron gas. Phys. Rev. B 1998, 58, 13465. (49) Wang, Y. A.; Govind, N.; Carter, E. A. Orbital-free kinetic-energy density functionals with a density-dependent kernel. Phys. Rev. B 1999, 60, 16350. (50) Shin, I.; Carter, E. A. Enhanced von Weizsäcker Wang-Govind-Carter kinetic energy density functional for semiconductors. J. Chem. Phys. 2014, 140, 18A531. (51) Ho, G. S.; Lignères, V. L.; Carter, E. A. Analytic form for a nonlocal kinetic energy functional with a density-dependent kernel for orbital-free density functional theory under periodic and Dirichlet boundary conditions. Phys. Rev. B 2008, 78, 045105. (52) Alonso, J. A.; Girifalco, L. A. Nonlocal approximation to the exchange potential and kinetic energy of an inhomogeneous electron gas. Phys. Rev. B 1978, 17, 3735. (53) Garc´ıa-González, P.; Alvarellos, J. E.; Chacón, E. Nonlocal kinetic-energy-density functionals. Phys. Rev. B 1996, 53, 9509.

34


Page 34 of 41

Page 35 of 41 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


(54) Chacón, E.; Alvarellos, J. E.; Tarazona, P. Nonlocal kinetic energy functional for nonhomogeneous electron systems. Phys. Rev. B 1985, 32, 7868. (55) Garc´ıa-González, P.; Alvarellos, J. E.; Chacón, E. Nonlocal symmetrized kineticenergy density functional: Application to simple surfaces. Phys. Rev. B 1998, 57, 4857. (56) Garcia-Aldea, D.; Alvarellos, J. E. Approach to kinetic energy density functionals: Nonlocal terms with the structure of the von Weizsäcker functional. Phys. Rev. A 2008, 77, 022502. (57) Lude˜ na, E. V.; Salazar, E. X.; Cornejo, M. H.; Arroyo, D. E.; Karasiev, V. V. The Liu-Parr power series expansion of the Pauli kinetic energy functional with the incorporation of shell-inducing traits: Atoms. Int. J. Quantum Chem. 2018, (58) Salazar, E. X.; Guarderas, P. F.; Ludena, E. V.; Cornejo, M. H.; Karasiev, V. V. Study of some simple approximations to the non-interacting kinetic energy functional. Int. J. Quantum Chem. 2016, 116, 1313–1321. (59) Perrot, F. Hydrogen-hydrogen interaction in an electron gas. J. Phys.: Cond. Mat. 1994, 6, 431. (60) Wang, L.-W.; Teter, M. P. Kinetic-energy functional of the electron density. Phys. Rev. B 1992, 45, 13196. (61) Smargiassi, E.; Madden, P. A. Orbital-free kinetic-energy functionals for firstprinciples molecular dynamics. Phys. Rev. B 1994, 49, 5220. (62) Garcia-Aldea, D.; Alvarellos, J. Fully nonlocal kinetic energy density functionals: A proposal and a general assessment for atomic systems. J. Chem. Phys. 2008, 129, 074103.

35



(63) Constantin, L. A.; Fabiano, E.; Laricchia, S.; Della Sala, F. Semiclassical neutral atom as a reference system in density functional theory. Phys. Rev. Lett. 2011, 106, 186406. (64) Laricchia, S.; Fabiano, E.; Constantin, L.; Della Sala, F. Generalized gradient approximations of the noninteracting kinetic energy from the semiclassical atom theory: Rationalization of the accuracy of the frozen density embedding theory for nonbonded interactions. J. Chem. Theory Comput. 2011, 7, 2439–2451. (65) Ou-Yang, H.; Levy, M. Approximate noninteracting kinetic energy functionals from a nonuniform scaling requirement. Int. J. Quantum Chem. 1991, 40, 379–388. (66) Perdew, J. P. Generalized gradient approximation for the fermion kinetic energy as a functional of the density. Phys. Lett. A 1992, 165, 79–82. (67) Ernzerhof, M. The role of the kinetic energy density in approximations to the exchange energy. J. Mol. Str.: THEOCHEM 2000, 501, 59–64. (68) Constantin, L. A.; Ruzsinszky, A. Kinetic energy density functionals from the Airy gas with an application to the atomization kinetic energies of molecules. Phys. Rev. B 2009, 79, 115117. (69) Vitos, L.; Johansson, B.; Kollar, J.; Skriver, H. L. Local kinetic-energy density of the Airy gas. Phys. Rev. A 2000, 61, 052511. (70) Lindmaa, A.; Mattsson, A. E.; Armiento, R. Quantum oscillations in the kinetic energy density: Gradient corrections from the Airy gas. Phys. Rev. B 2014, 90, 075139. ´ (71) Constantin, L. A.; Fabiano, E.; Smiga, S.; Della Sala, F. Jellium-with-gap model applied to semilocal kinetic functionals. Phys. Rev. B 2017, 95, 115153. (72) Borgoo, A.; Tozer, D. J. Density scaling of noninteracting kinetic energy functionals. J. Chem. Theory Comput. 2013, 9, 2250–2255.

36


Page 36 of 41

Page 37 of 41 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


(73) Thakkar, A. J. Comparison of kinetic-energy density functionals. Phys. Rev. A 1992, 46, 6920. (74) Lembarki, A.; Chermette, H. Obtaining a gradient–corrected kinetic–energy functional from the Perdew–Wang exchange functional. Phys. Rev. A 1994, 50, 5328. (75) Xia, J.; Carter, E. A. Single-point kinetic energy density functionals: A pointwise kinetic energy density analysis and numerical convergence investigation. Phys. Rev. B 2015, 91, 045124. (76) Thomas, L. H. The calculation of atomic fields. Math. Proc. Cambridge Phil. Soc. 1927; pp 542–548. (77) Fermi, E. Un metodo statistico per la determinazione di alcune priorieta dellatome. Rend. Accad. Naz. Lincei 1927, 6, 32. (78) von Weizsäcker, C. F. Zur theorie der kernmassen. Zeitschrift f¨ ur Physik A Hadrons and Nuclei 1935, 96, 431–458. (79) Levy, M.; Ou-Yang, H. Exact properties of the Pauli potential for the square root of the electron density and the kinetic energy functional. Phys. Rev. A 1988, 38, 625. (80) Lude˜ na, E. V.; Arroyo, D.; Salazar, E. X.; Vallejo, J. In Novel Electronic Structure Theory: General Innovations and Strongly Correlated Systems; Hoggan, P. E., Ed.; Advances in Quantum Chemistry; Academic Press, 2018; Vol. 76; pp 59 – 78. (81) Kirzhnitz, D. Quantum corrections to the Thomas–Fermi equation. Sov. Phys. JETP 1957, 5, 64. (82) Yan, Z.; Perdew, J. P.; Korhonen, T.; Ziesche, P. Numerical test of the sixth-order gradient expansion for the kinetic energy: Application to the monovacancy in jellium. Phys. Rev. A 1997, 55, 4601.

37



(83) Hodges, C. H. Quantum corrections to the Thomas–Fermi approximation: the Kirzhnits method. Can. J. Phys. 1973, 51, 1428. (84) Brack, M.; Jennings, B. K.; Chu, Y. H. On the extended Thomas–Fermi approximation to the kinetic energy density. Phys. Lett. B 1976, 65, 1–4. (85) Constantin, L. A.; Fabiano, E.; Della Sala, F. Modified fourth-order kinetic energy gradient expansion with an Hartree potential dependent correction. J. Chem. Theory Comput. 2017, 13, 4228–4239. (86) Zhou, B.; Alexander Wang, Y.; Carter, E. A. Transferable local pseudopotentials derived via inversion of the Kohn-Sham equations in a bulk environment. Phys. Rev. B 2004, 69, 125109. (87) Huang, C.; Carter, E. A. Transferable local pseudopotentials for magnesium, aluminum and silicon. Phys. Chem. Chem. Phys. 2008, 10, 7109–7120. (88) Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 1981, 23, 5048. (89) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (90) Jana, S.; Patra, A.; Samal, P. Assessing the performance of the Tao-Mo semilocal density functional in the projector-augmented-wave method. J. Chem. Phys. 2018, 149, 044120. (91) Xia, J.; Carter, E. A. Density-decomposed orbital-free density functional theory for covalently bonded molecules and materials. Phys. Rev. B 2012, 86, 235109. (92) Laricchia, S.; Fabiano, E.; Della Sala, F. On the accuracy of frozen density embedding calculations with hybrid and orbital-dependent functionals for non-bonded interaction energies. J. Chem. Phys. 2012, 137, 014102. 38


Page 38 of 41

Page 39 of 41 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


(93) Xia, J.; Huang, C.; Shin, I.; Carter, E. A. Can orbital-free density functional theory simulate molecules? J. Chem. Phys. 2012, 136, 084102. (94) Iyengar, S. S.; Ernzerhof, M.; Maximoff, S. N.; Scuseria, G. E. Challenge of creating accurate and effective kinetic-energy functionals. Phys. Rev. A 2001, 63, 052508. (95) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (96) Gonze, X. First-principles computation of material properties: the ABINIT software project. Comp. Mater. Sci. 2002, 25, 478–492. (97) Birch, F. Finite Elastic Strain of Cubic Crystals. Phys. Rev. 1947, 71, 809–824. (98) Elliott, P.; Lee, D.; Cangi, A.; Burke, K. Semiclassical origins of density functionals. Phys. Rev. Lett. 2008, 100, 256406. (99) Lee, D.; Constantin, L. A.; Perdew, J. P.; Burke, K. Condition on the Kohn–Sham kinetic energy and modern parametrization of the Thomas–Fermi density. J. Chem. Phys. 2009, 130, 034107. (100) Elliott, P.; Burke, K. Non-empirical derivation of the parameter in the B88 exchange functional. Can. J. Chem. 2009, 87, 1485–1491. (101) Constantin, L. A.; Terentjevs, A.; Della Sala, F.; Cortona, P.; Fabiano, E. Semiclassical atom theory applied to solid-state physics. Phys. Rev. B 2016, 93, 045126. (102) Laricchia, S.; Fabiano, E.; Constantin, L. A.; Della Sala, F. Generalized Gradient Approximations of the Noninteracting Kinetic Energy from the Semiclassical Atom Theory: Rationalization of the Accuracy of the Frozen Density Embedding Theory for Nonbonded Interactions. J. Chem. Theory Comput. 2011, 7, 2439–2451. (103) Fabiano, E.; Constantin, L. A. Relevance of coordinate and particle-number scaling in density-functional theory. Phys. Rev. A 2013, 87, 012511. 39



(104) Constantin, L. A.; Snyder, J. C.; Perdew, J. P.; Burke, K. Communication: Ionization potentials in the limit of large atomic number. J. Chem. Phys. 2011, 133, 241103.

40


Page 40 of 41

Page 41 of 41 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


Graphical TOC Entry

41


Performance of semilocal kinetic-energy functionals for Orbital-Free

Recommend Documents