Chemical Physics Letters

(1)

Orbital-free density functional theory: Linear scaling methods for kinetic potentials, and applications to solid Al and Si

Jeng-Da Chai

^a,1

, Vincent L. Lignères

^b

, Gregory Ho

^b

, Emily A. Carter

^c

, John D. Weeks

^a,d,*

aInstitute for Physical Science and Technology, University of Maryland, College Park, MD 20742, United States

bDepartment of Chemistry, Princeton University, Princeton, NJ 08544, United States

cDepartment of Mechanical and Aerospace Engineering and Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, United States

dDepartment of Chemistry and Biochemistry, University of Maryland, College Park, MD 20742, United States

a r t i c l e i n f o

Article history:

Received 11 October 2008 In ﬁnal form 23 March 2009 Available online 27 March 2009

a b s t r a c t

In orbital-free density functional theory the kinetic potential (KP), the functional derivative of the kinetic energy density functional, appears in the Euler equation for the electron density and may be more ame- nable to simple approximations. We study properties of two solid-state systems, Al and Si, using two nonlocal KPs that gave good results for atoms. Very accurate results are found for Al, but results for Si are much less satisfactory, illustrating the general need for a better treatment of extended covalent systems. A different integration pathway in the KP formalism may prove useful in attacking this fundamen- tal problem.

1. Introduction

Density-functional theory (DFT) is one of the most useful methods for studying ground state properties of many-electron systems.

In principle, as shown by Hohenberg and Kohn [1], the exact ground state energy of a system of N electrons can be formally written as a functional E½

q

of only the electron density

q

ðrÞ, a function of three variables, and the external ﬁeld VextðrÞ[2,3]. Kohn and Sham (KS)[4]partitioned E½

q

into the sum of three terms, the classical Hartree energy EH½

q

(readily expressed in terms of the density), the kinetic energy density functional (KEDF) Ts½

q

of a model system of noninteracting electrons at the same density, and the much smaller remainder, the so-called exchange–correlation energy functional Exc½

q

. KS showed that the numerical value of Ts½

q

could be determined exactly, not directly from the density itself, but by using a set of N one-electron wave functions (orbitals) to solve the N coupled KS equations that describe the noninteracting system. By using these results along with relatively simple approximations for the small Exc½

q

term, the resulting KS-DFT has proven successful in many applications.

However, because of the use of the KS orbitals, KS-DFT typically scales as OðN³Þ, which limits the applicability of KS-DFT. By using

localized orbitals some orbital-based linear-scaling methods have been developed that are especially successful for insulating systems with large band gaps[5,6]. However, the prefactors of these types of approaches are usually so large that their computational costs become cheaper than the traditional KS-DFT (with OðN³Þ scaling) only for very large systems.

This problem could be circumvented if we could develop an accurate ‘orbital-free’ density functional theory (OF-DFT), where the kinetic energy is expressed in terms of the electron density alone[2,3,7,8]. Of course, given its large magnitude, this is a very difﬁcult task, and indeed simple local approximations like that used in the Thomas–Fermi (TF) model [9,10] have proved very inaccurate.

Through the recent efforts of many workers[11–23], we now have significantly better approximations for the KEDF. There have been two main advances. The first is the use of nonlocal KEDFs that reproduce known exact results both for very slowly varying or very rapidly varying fields, and that give the exact linear response (LR) of the density of uniform model systems to small perturbations.

The second is to focus on the more slowly-varying valence electron density as described by weak local pseudopotentials[24–34]. For OF-DFT with LR-based KEDFs, the use of such pseudopotentials not only can reduce the computational cost relative to all-electron calculations, but also can improve its accuracy, since the valence system will be closer to the LR regime where T_s½

q

is designed to be accurate [24,33]. Very promising results using such OF-DFT methods have been obtained for a variety of nearly free-electron- like metals.

However, a more accurate treatment of the KEDF is still needed in other applications with signiﬁcant variations in the density. The 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.

doi:10.1016/j.cplett.2009.03.064

* Corresponding author. Address: Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, United States.

E-mail addresses: [email protected] (J.-D. Chai), [email protected] (J.D.

Weeks).

1 Present address: Department of Chemistry, University of California at Berkeley, and Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States.

Contents lists available atScienceDirect

Chemical Physics Letters

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c p l e t t

(2)

main problem is that the exact Ts½

q

is highly nonlocal, and we have little idea of the functional form of the nonlocality for densities far from the linear response regime.

In earlier work, Chai and Weeks (CW)[34,35]proposed a different way to attack this basic problem. In the usual approach, approximate forms for Ts½

q

and Exc½

q

are assumed, and the density

q

ðrÞ is obtained from the variational principle (Euler equation) associated with minimizing E½

q

:

l

¼ VTsðr; ½

q

Þ þ Veffðr; ½

q

Þ: ð1Þ

The total energy of the inhomogeneous system is then determined from the energy functional E½

q

. All other physical quantities related to the ground-state density could also be computed. Here

l

is the chemical potential (the Lagrange multiplier associated with the normalization condition R

q

ðrÞdr ¼ N), and Veffðr; ½

q

Þ is an effective one-body potential deﬁned by

Veffðr; ½

q

Þ VHðr; ½

q

Þ þ Vxcðr; ½

q

Þ þ VextðrÞ; ð2Þ where

VHðr; ½

q

Þ Z

q

ðr⁰Þ

jr r⁰jdr⁰ ð3Þ

is the Hartree potential, and Vxcðr; ½

q

Þ dExc½

q

=d

q

ðrÞ is the exchange–correlation potential. Similarly we interpret

VTsðr; ½

q

Þ dTs½

q

=d

q

ðrÞ ð4Þ

as the kinetic potential (KP) arising from the KEDF[36,37].

In principle this process could be reversed, and the exact Ts½

q

could be determined from the KP VTsðr; ½

q

Þ by functional integration over density changes in all regions of space[22,34,36,37]. Be- cause of this integration, Ts½

q

is a functional that depends on the density everywhere even in the simplest case, exempliﬁed by TF theory, where VTsðr; ½

q

Þ is a simple function of the local density at r only.

Of course, the true VTsðr; ½

q

Þ itself is much more complicated and must be a nonlocal functional of the density as well. However it seems reasonable to assume that it depends most strongly on the density in some local region around r. Recent detailed arguments [38,39]suggest that the KP in this sense is more local than is the KEDF. Since most problems in devising accurate approximations for Ts½

q

have arisen from the nonlocality, this suggests that it could be worthwhile to try to develop approximations for the more local KP VTsðr; ½

q

Þ.

To that end, following ideas ﬁrst introduced for the OF-KEDF’s, CW proposed two different nonlocal KPs for atomic systems that reduce to known exact forms for both slowly varying and rapidly varying perturbations and also reproduce exact results for the linear response of the density of the homogeneous system to small perturbations[34]. The latter requirement introduces nonlocality into the resulting KPs. CW also proposed two efﬁcient pathways to determine the numerical value of the kinetic energy by integration of the KPs, as discussed below in Section3.

Using the CW KPs and the most accurate (Herring) pathway along with weak ab initio local pseudopotentials for the valence component of the electron density, CW obtained results for the total energies and ionization energies of atoms, and for the shell structure in the atomic radial density proﬁles that are in better agreement with calculations using the full Kohn–Sham theory than earlier KEDFs were able to produce. For example, the ionization energies of the ﬁrst and second row atoms predicted by the CW KPs have average errors of 1.8 eV and 2.1 eV for the two different KPs, compared with the full KS-DFT, while the CAT KEDF [17], a nonlocal KEDF, has an average error of 3.9 eV. Similarly, errors in the valence energies predicted by CAT KEDF are almost twice as large than those predicted by the CW KPs[34].

In this Letter, we use the CW KPs to study two very different solid-state systems, a nearly-free-electron-like metal (solid Al) and a covalent material (solid Si). As discussed below, the original forms of these KPs use the local Fermi wave vector. This permits their use in atomic and molecular systems where the electron density van- ishes far from the nuclei. However, this also leads to a quadratic scaling with system size[34], which limits their usefulness for very large extended systems. Following earlier work on the KEDF [12,14], we introduce in Section2an expansion method to reduce the computational cost of evaluating these nonlocal KPs, and the truncated expansions essentially scale linearly with system size.

The density pathway by which the kinetic energy is obtained from the KP is described in Section3. Section4compares the results of this method with the full KS-DFT and with other approaches using KEDFs and gives our conclusions.

2. Linear scaling methods

CW introduced two different nonlocal KPs [34], one denoted LQ because it satisﬁes known exact results to second order at low wavevectors (‘Low q’), and the other HQ with similar behav- ior at high wavevectors (‘High q’). The nonlocality is generated by the requirement that these KPs also reproduce the exact linear response function of a uniform electron gas at the local Fer- mi wave vector (LFWV) kFðrÞ, deﬁned in terms of the electron density by

kFðrÞ ð3

p

²

q

ðrÞÞ¹⁼³: ð5Þ

Both potentials can be compactly described using the following generalized form with different values for a parameter

a

: Vaðr; ½

q

; kFðrÞÞ ¼ VTFðr; ½

q

Þ þ VWðr; ½

q

Þ þ10

9

a

^C^F

q

^bðrÞ

Z

f ðjr r⁰j; kFðrÞÞ

q

^aðr⁰Þdr⁰ ð6Þ Here b ð2=3

a

Þ. The HQ model uses

a

¼ 2=3 and the LQ model has

a

¼ 1=2. VTFðr; ½

q

Þ and VWðr; ½

q

Þ are the TF KP and the von Weizsäcker (W) KP[40]respectively, and C_Fis a numerical coefﬁ- cient associated with the TF KEDF. See Ref. [34] for a detailed discussion.

The weight function f ðjr r⁰j; kFðrÞÞ is deﬁned by the following inverse Fourier-transform-like integral:

f ðjr r⁰j; kFðrÞÞ 1 ð2

p

Þ³

Z ^f ðk=2kFðrÞÞe^ikðrr⁰^Þdk: ð7Þ

Here

^f ðqÞ ¼ FLðqÞ 3q² 1: ð8Þ

is directly related to the inverse linear response function FLðqÞ of the uniform electron gas with density

q

₀ at reduced wave vector q k=2kF0, with kF0 ð3

p

²

q

₀Þ¹⁼³[41]:

FLðqÞ 1 2þ1 q²

4q ln 1 þ q 1 q

¹

: ð9Þ

Clearly kF0is replaced by the LFWV in the the deﬁnition of the weight function f ðjr r⁰j; kFðrÞÞ in Eq.(7). Given this deﬁnition it is straightforward to rewrite the last term in Eq. (6) in Fourier space as

10

9

a

ð2

p

Þ³CF

q

^bðrÞ

Z ^f ðk=2kFðrÞÞ

q

^aðkÞe^ikrdk: ð10Þ

Because of the kFðrÞ term in Eq.(10), the integration must be done on a grid over a range of r, and all these models scale quadratically.

(3)

In order to achieve the desirable linear scaling for fast computa- tion in large solid-state systems, we follow the work of Wang et al.

[12]and apply a Taylor series expansion to the weight functional f ðjr r⁰j; kFðrÞÞ with respect to a reference density

q

. For extended systems, where

q

ðrÞ is not signiﬁcantly different from the average density

q

₀, the natural choice of

q

is

q

₀. Throughout this Letter, we choose

q

¼

q

₀in all of our calculations.

The weight function f ðjr r⁰j; kFðrÞÞ can be expanded with respect to

q

as:

f ðjr r⁰j; kFðrÞÞ ¼ f ðjr r⁰j; kFÞ þ f^ð1Þðjr r⁰j; kFÞð

q

ðrÞ

q

Þ þ1

2f^ð2Þðjr r⁰j; kFÞð

q

ðrÞ

q

Þ²þ ð11Þ where k_F ð3

p

²

q

Þ¹⁼³, and f^ðnÞðjrr⁰j;k_FÞ¼@ⁿf ðjrr⁰j;kFðrÞÞ=@

q

ⁿðrÞj_q is the nth order derivative of f ðjr r⁰j; kFðrÞÞ with respect to

q

ðrÞ, and is evaluated at

q

. Their functional forms, up to second order, in reciprocal space are

^f^ð1ÞðqÞ ¼ q 3

q

^{^}^f

0ðqÞ ð12Þ

and

^f^ð2ÞðqÞ ¼q²^f⁰⁰ðqÞ þ 4q^f⁰ðqÞ

ð3

q

Þ² ð13Þ

where q¼ k=ð2k_FÞ, and ^f⁰ðqÞ and ^f⁰⁰ðqÞ are the ﬁrst and the second derivative of ^fðqÞ with respect to q. Since the analytical form of ^fðqÞ is available in Eq. (8), all the terms needed in the Taylor series expansion ^f^ðnÞðqÞ can be obtained analytically. This simplicity is one of the main advantages of the KP method over related methods that use the KEDF, where these terms have usually been be obtained numerically, with one recent exception[42].

By carrying out the Taylor series expansion to nth order in Eq.

(11), and inserting it in Eq. (6), one can then take out the ð

q

ðrÞ

q

Þ^mfactor of the mth term from its integral (where m is a nonnegative integer, and m 6 n), and the remaining integral becomes:

Z

f^ðmÞðjr r⁰j; kFÞ

q

^aðr⁰Þdr⁰ ð14Þ The integral in Eq.(14)can be easily evaluated by a simple fast Fourier transform (FFT). Therefore, to compute the nonlocal term in Eq. (6) using this scheme, one needs to evaluate a total of n þ 2 FFT’s, including the FFT of

q

^a. Since all of the terms can be computed by FFT’s, this scheme essentially scales linearly OðM ln MÞ with system size, where M is the number of grid points.

The corresponding linear-scaling HQ ð

a

¼ 2=3Þ and LQ ð

a

¼ 1=2Þ KPs with

q

¼

q

₀can be written as:

V^lin_aðr; ½

q

; kF0Þ ¼ VTFðr; ½

q

Þ þ VWðr; ½

q

Þ þ10

9

a

^C^F

q

^bðrÞ Z

f ðjr r⁰j; kF0Þ

q

^aðr⁰Þdr⁰þ ð

q

ðrÞ

q

₀Þ

Z

f^ð1Þðjr r⁰j; kF0Þ

q

^aðr⁰Þdr⁰þ1

2ð

q

ðrÞ

q

0Þ²

Z

f^ð2Þðjr r⁰j; kF0Þ

q

^aðr⁰Þdr⁰þ

ð15Þ where kF0 ð3

p

²

q

₀Þ¹⁼³is the uniform Fermi wave vector. For simple nearly-free-electron metals like Al, the expansion method works very well, and indeed often only the zeroth order term is needed.

For extended systems with large density variations over space, this expansion method could experience convergence problems similar to those that have been seen by other workers[12–14,24,25].

However, in those cases, both the use of the LFWV in Eq.(5)and the basic linear response treatment of the nonlocality are probably inadequate as well. A two-body Fermi wavevector kFðr; r⁰Þ like the one used in the WGC KEDF[12]could be introduced, but this func-

tional form represents an additional approximation that does not systematically improve the underlying linear response treatment of the nonlocality. As we will see, the treatment of extended systems with large density variations remains a major challenge for all OF-DFT methods.

3. Density pathway

As discussed in Ref.[34], an integration pathway is needed to determine the value of Ts½

q

from a given VTsðr; ½

q

Þ [43–46]. If the exact VTsðr; ½

q

Þ is used and the integration is carried out exactly, then all pathways would give the same exact result for Ts½

q

. An approximate KP can give the same (approximate) value for the kinetic energy independent of pathway only if it arises from functionally differentiating a KEDF (the usual OFDFT approach) or equivalently if it exactly satisﬁes the nonlocal consistency conditions given by Herring[22]in his Eq. (24). As he points out, it is not easy to satisfy these formal conditions for general functionals, and we do not try to do so here.

Instead, we generate nonlocality by imposing a more easily implemented and physically suggestive condition, requiring that linear response theory is exactly satisﬁed for small perturbations about a uniform system. Thus our results will be independent of pathway in the linear regime but the kinetic energy obtained from integration of our approximate KPs for systems with large density variations will in general depend on the particular integration pathway used, as discussed in detail in Ref.[34]. Experience with the classical analogue of the KP method applied to nonuniform hard sphere ﬂuids [46] has shown that errors from this path dependence can be small (e.g., less than one percent for the surface tension at a hard wall) when used with particular pathways that do not excessively weight regions poorly described by an approximate theory for the nonuniform density.

Two efﬁcient pathways to obtain the kinetic energy from an approximate KP for atomic and molecular systems were discussed in Ref.[34]. The simplest and most accurate pathway, due to Her- ring [22], automatically satisﬁes the virial theorem. However it does not apply to extended systems in its original form and we have not yet found the appropriate generalization. Here we use the alternate linear density pathway, where there is a linear scaling of the electron density

q

_kðrÞ

q

₀þ k½

q

ðrÞ

q

₀: ð16Þ

If the exact KP were used, both pathways would give the exact kinetic energy.

By inserting the HQ and LQ KPs of Eq.(15)into the linear density pathway (see Eqs.(9)–(11)in Ref.[34]), the kinetic energy can be computed as

Ta¼ Tk¼0þ Z 1

0

dk Z

dr½

q

ðrÞ

q

₀V^lin_aðr; ½

q

_k; kF0Þ: ð17Þ Here T_k¼0is the kinetic energy for uniform system, i.e., the Thomas–

Fermi kinetic energy TTF½

q

₀. Since the VTFðr; ½

q

_kÞ and VWðr; ½

q

_kÞ in Eq.(17)arise from the functional derivatives of the known TTF½

q

and TW½

q

functionals respectively, these terms can be integrated exactly and lead to the TF and W KEDFs.

Eq.(17)then becomes Ta½

q

¼ TTF½

q

þ TW½

q

þ

Z 1 0

dk Z

dr½

q

ðrÞ

q

₀

10 9

a

^C^F

q

^b_kðrÞ

Z

dr⁰f ðjr r⁰j; kF0Þ

q

^a_kðr⁰Þ þ ð

q

_kðrÞ

q

₀Þ

Z

dr⁰f^ð1Þðjr r⁰j; kF0Þ

q

^a_kðr⁰Þ þ1

2ð

q

_kðrÞ

q

₀Þ²

Z

dr⁰f^ð2Þðjr r⁰j; kF0Þ

q

^a_kðr⁰Þ þ

ð18Þ

(4)

where Tk¼0¼ TTF½

q

₀ term is absorbed in the TTF½

q

term. We use Eq.

(18) to evaluate the kinetic energy of our linear-scaling HQ ð

a

¼ 2=3Þ and LQ ð

a

¼ 1=2Þ KPs in all the calculations here.

4. Results

In solids, the external potential VextðrÞ can be regarded as a linear combination of the special array of local atomic pseudopotentials centered at each ion position RI. Different arrays of RIlead to different phases, such as face-centered cubic (fcc), diamond (dia), body-centered cubic (bcc), simple cubic (sc), and so on. We calcu- late the binding energies of Al and Si at these four different phases, and compare our results with other KEDFs and KS-DFT. We use periodic boundary conditions with a cubic supercell containing 4 atoms for fcc, 8 atoms for dia, 2 atoms for bcc, and 8 atoms for sc. All calculations are spin-restricted and use the local density approximation (LDA) [47–50] for the exchange–correlation functional.

For bulk Al, the empirical Goodwin–Needs–Heine (GNH) local pseudopotential[29]is used, and a plane wave kinetic energy cutoff of 600 eV is used to converge the electron density. Here it is suf- ficient to use only the zeroth-order linear-scaling HQ and LQ KPs (see Eq.(15)) in the calculations, as verified by comparison with results of a first order calculation. Our results are compared with the WT[23]and WGC[12]KEDFs, and KS-DFT, which were previously computed[12], using the same local pseudopotential. (The atomic- based local pseudopotentials in Ref. [34]could also be used and gives very similar results.) As can be seen inTables 1 and 2, all the LR-based models perform similarly, and agree well with KS- DFT. The phase ordering is correct, and aside from the high energy diamond phase, the lattice parameters are close to the KS results.

To assess the performance of the HQ and LQ KPs in covalent systems like solid Si, we use the bulk local pseudopotential (BLPS) developed by Zhou et al.[24]together with a plane wave kinetic energy cutoff of 2000 eV for converging the electron density. We used the ﬁrst-order linear-scaling HQ and LQ KPs in the Si calculations. Our results are compared ﬁrst with the original WGC KEDF [12], and with the KS-DFT, which were previously computed [24], using the same local pseudopotential. As can be seen inTables 3 and 4, all these LR-based models perform worse than for Al, when compared with KS-DFT. The phase ordering is incorrect, and the lattice parameters of the four phases only qualitatively match with the KS results, although the WGC KEDF at least obtains the correct diamond ground state.

Because we use approximate KPs, our estimate for Ts½

q

will depend on the integration pathway, and it is important to try to find a particular pathway that is relatively insensitive to the errors that exist in our KPs. For atomic systems, we were able to use the Her- ring pathway, which automatically satisfies the virial theorem, and good results were found[34]. As noted above, a generalization of the Herring pathway for solids is not yet available, and in this Let- ter we used the linear density pathway. While this pathway is numerically efficient, and gives exact results when the exact KP is used, we expect less accurate results when used with our approximate KPs since the virial theorem is not satisfied.

Indeed this is the case for atomic systems, where results from both pathways can be compared. InTable 5, we show the pathway dependence of the full LQ KP (with ab initio local pseudopotentials [34]) in atomic systems. The kinetic energies of the three atoms (Al, Si and Ar), using the full LQ KP, are evaluated by the Herring and the density pathways, and are compared with the KS results.

Clearly, the density pathway gives signiﬁcantly worse kinetic energy than the Herring’s pathway, and its errors increase when the densities of the systems are more rapidly varying. This suggests that the large deviations of the LQ and HQ results for solid Si (dia) inTable 4may be partially due to the use of this pathway. This pathway dependence also seems likely to inﬂuence the phase orderings. Thus the generalization of the accurate Herring pathway to solids seems a promising focus for future research.

Since Si is a covalent material, the density inside the covalent bond regions is quite different from that outside. Even with a local pseudopotential this system may be outside the linear response regime, or the results at least may depend sensitively of the choice of the one- or two-body FWV. In particular, recent work [14] has shown that much improved results for the phase energies can be obtained for bulk Si using the WGC KEDF with a FWV mixing Table 1

Lattice parameters (Å) for bulk Al. The results for the KS, WT, and WGC models are taken from Ref.[12].

Al KS LQ HQ WT WGC

fcc 4.03 4.04 4.04 4.04 4.03

bcc 3.23 3.23 3.23 3.23 3.22

sc 5.33 5.31 5.36 5.33 5.38

dia 5.84 5.91 5.90 5.94 5.92

Table 2

Energy per atom (eV) for bulk Al. The ﬁrst row is the energy for the fcc structure, while other rows are energy difference from the fcc structure. The results for the KS, WT, and WGC models are taken from Ref.[12].

Al KS LQ HQ WT WGC

fcc 58.336 58.303 58.314 58.331 58.331

bcc 0.068 0.053 0.057 0.060 0.066

sc 0.250 0.253 0.253 0.227 0.217

dia 0.599 0.712 0.751 0.673 0.584

Table 3

Lattice parameters (Å) for bulk Si. The results for the KS and WGC models are taken from Ref.[24], and those for the WGC2 model are taken from Ref.[14].

Si KS LQ HQ WGC WGC2

dia 5.38 5.27 5.29 5.77 5.57

bcc 3.29 3.06 3.06 3.29 3.32

sc 4.99 4.98 4.98 5.01 5.06

fcc 3.83 3.83 3.82 3.80 3.80

Table 4

Energy per atom (eV) for bulk Si. The ﬁrst row is the energy for the dia structure, while other rows are energy difference from the dia structure. The results for the KS and WGC models are taken from Ref.[24], and those for the WGC2 model are taken from Ref.[14].

Si KS LQ HQ WGC WGC2

dia 110.234 109.167 109.282 110.345 110.220

bcc 0.165 0.553 0.437 0.537 0.267

sc 0.303 0.586 0.531 0.506 0.308

fcc 0.457 0.584 0.478 0.571 0.444

Table 5

The kinetic energy Ts½q (eV) of the three atoms (Al, Si, and Ar) in ab initio local pseudopotential calculations[34], using the KS method and the full LQ model (evaluated by both the Herring’s and the density pathways). MAE, the mean absolute errors (relative to the KS method) of the full LQ model computed by the two pathways are given at the bottom of their respective columns.

KS LQH LQden

Al 21.164 21.823 17.563

Si 40.985 41.264 33.874

Ar 225.609 221.035 196.495

MAE 1.838 13.276

(5)

parameter and reference density optimized for the covalent phases of Si, rather than the original FWV mixing parameter optimized for Al[12]. Results with this new parameter set, denoted here WGC2, are also given inTables 3 and 4, where we see that the phase orderings and energy differences for Si phases using WGC2 are in very good agreement with KS-DFT. However, ultimately it is desirable to ﬁnd a general KEDF or KP that does not depend on special properties of the system, particularly when considering more complicated systems with defects or surfaces. In our opinion, developing a truly universal and yet accurate KEDF or KP for covalent materials remains an outstanding problem.

5. Conclusion

In summary, CW have previously demonstrated that the nonlocal HQ and LQ KPs work well for isolated atoms and ions[34], and now we show they can be used with no additional parame- terization in solid-state systems. Very good results are found for Al. However, problems of inaccurate phase ordering and bond lengths in Si are found, and a better treatment of nonlocality be- yond the linear response regime may be needed for covalent systems with signiﬁcant density variations. The simpler and more local form of the OF KP could be useful in future developments.

If it is possible to develop a new pathway to obtain the kinetic energy in extended systems that satisﬁes the virial theorem, anal- ogous to the accurate Herring pathway for atomic systems[34], we believe even the present LR-based KP’s would likely give more accurate phase energies. Further work along these lines is called for.

Acknowledgements

This work at the University of Maryland has been supported by the NSF Grants No. CHE01-11104 and CHE05-17818, and by the NSF-MRSEC under Grants No. DMR 00-80008 and DMR05-20471.

E.A.C. is grateful for support from the NSF under Grant No.

CHE05-17359 for this work. One of the authors (J.D.C.) acknowl- edges the support from the UMCPGraduate School program, the IPST

Alexander program, and the CHPH Block Grant Supplemental program.

References

[1] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864.

[2] R.G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, 1989.

[3] R.M. Dreizler, E.K.U. Gross, Density Functional Theory: An Approach to the Quantum Many Body Problem, Springer-Verlag, Berlin, 1990.

[4] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133.

[5] W. Yang, Phys. Rev. Lett. 66 (1991) 1438.

[6] For a recent review, see S. Goedecker, Rev. Mod. Phys. 71 (1999) 1085.

[7] S.C. Watson, E.A. Carter, Comput. Phys. Commun. 128 (2000) 67.

[8] N. Choly, E. Kaxiras, Solid State Commun. 121 (2002) 281.

[9] L.H. Thomas, Proc. Cambridge Phil. Soc. 23 (1927) 542.

[10] E. Fermi, Z. Phys. 48 (1928) 73.

[11] Y.A. Wang, N. Govind, E.A. Carter, Phys. Rev. B 58 (1998) 13465;

Y.A. Wang, N. Govind, E.A. Carter, Phys. Rev. B 60 (1999) 17162(E);

Y.A. Wang, N. Govind, E.A. Carter, Phys. Rev. B 64 (2001) 129901(E).

[12] Y.A. Wang, N. Govind, E.A. Carter, Phys. Rev. B 60 (1999) 16350;

Y.A. Wang, N. Govind, E.A. Carter, Phys. Rev. B 64 (2001) 089903(E).

[13] See e.g. Y.A. Wang, E.A. Carter, in: S.D. Schwartz (Ed.), Theoretical Methods in Condensed Phase Chemistry, Progress in Theoretical Chemistry and Physics, Kluwer, Boston, 2000, p. 117, and references therein.

[14] B. Zhou, V.L. Ligneres, E.A. Carter, J. Chem. Phys. 122 (2005) 044103.

[15] E. Chacón, J.E. Alvarellos, P. Tarazona, Phys. Rev. B 32 (1985) 7868.

[16] P. García-González, J.E. Alvarellos, E. Chacón, Phys. Rev. A 54 (1996) 1897.

[17] P. García-González, J.E. Alvarellos, E. Chacón, Phys. Rev. B 57 (1998) 4857.

[18] M. Foley, P.A. Madden, Phys. Rev. B 53 (1996) 10589.

[19] F. Perrot, J. Phys.: Condens. Matter 6 (1994) 431.

[20] J.A. Alonso, L.A. Girifalco, Phys. Rev. B 17 (1978) 3735.

[21] M.D. Glossman, L.C. Balbás, J.A. Alonso, Chem. Phys. 196 (1995) 455.

[22] C. Herring, Phys. Rev. A 34 (1986) 2614.

[23] L.-W. Wang, M.P. Teter, Phys. Rev. B 45 (1992) 13196.

[24] B. Zhou, Y.A. Wang, E.A. Carter, Phys. Rev. B 69 (2004) 125109.

[25] B. Zhou, E.A. Carter, J. Chem. Phys. 122 (2005) 184108.

[26] W.C. Topp, J.J. Hopﬁeld, Phys. Rev. B 7 (1973) 1295.

[27] J.A. Appelbaum, D.R. Hamann, Phys. Rev. 8 (1973) 1777.

[28] M. Schlüter, J.R. Chelikowsky, S.G. Louie, M.L. Cohen, Phys. Rev. 12 (1975) 4200.

[29] L. Goodwin, R.J. Needs, V. Heine, J. Phys.: Condes. Matter 2 (1990) 351.

[30] S. Watson, B.J. Jesson, E.A. Carter, P.A. Madden, Europhys. Lett. 41 (1998) 37.

[31] J.A. Anta, P.A. Madden, J. Phys.: Condens. Matter 11 (1999) 6099.

[32] D.J. González, L.E. González, J.M. López, M.J. Stott, Phys. Rev. B 65 (2002) 184201.

[33] B. Wang, M.J. Stott, Phys. Rev. B 68 (2003) 195102.

[34] J.-D. Chai, J.D. Weeks, Phys. Rev. B 75 (2007) 205122.

[35] J.-D. Chai, J.D. Weeks, J. Phys. Chem. B 108 (2004) 6870.

[36] R.A. King, N.C. Handy, Phys. Chem. Chem. Phys. 2 (2000) 5049.

[37] R.A. King, N.C. Handy, Mol. Phys. 99 (2001) 1005.

[38] A. Holas, N.H. March, Phys. Rev. A 66 (2002) 066501.

[39] I. Lindgren, S. Salomonson, Phys. Rev. A 67 (2003) 056501.

[40] C.F. von Weizsäcker, Z. Physik 96 (1935) 431.

[41] J. Lindhard, K. Dan Vidensk Selsk, Mat. Fys. Medd. 28 (1954) 8.

[42] G. Ho, V.L. Ligneres, E.A. Carter, Phys. Rev. B 78 (2008) 045105.

[43] L.R. Pratt, G.G. Hoffman, R.A. Harris, J. Chem. Phys. 88 (1988) 1818.

[44] L.R. Pratt, G.G. Hoffman, R.A. Harris, J. Chem. Phys. 92 (1990) 6687.

[45] G.G. Hoffman, L.R. Pratt, Mol. Phys. 82 (1994) 245.

[46] Y.-G. Chen, J.D. Weeks, J. Chem. Phys. 118 (2003) 7944.

[47] P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26 (1930) 376.

[48] D.M. Ceperley, Phys. Rev. B 18 (1978) 3126.

[49] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566.

[50] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048.