Nonlocal density scheme for electronic-structure calculations

1 downloads 0 Views 81KB Size Report
ids, while the LDA became the state-of-the-art of total- energy and band-structure calculations. Keeping into account that several realistic parametrizations of Kxc.
PHYSICAL REVIEW B

VOLUME 60, NUMBER 16

15 OCTOBER 1999-II

Nonlocal density scheme for electronic-structure calculations Maurizia Palummo, Giovanni Onida, Rodolfo Del Sole, and Massimiliano Corradini Istituto Nazionale per la Fisica della Materia–Dipartimento di Fisica dell’ Universita´ di Roma Tor Vergata, Via della Ricerca Scientifica, I-00133 Roma, Italy

Lucia Reining

Laboratoire des Solides Irradie´s, UMR7642 CNRS-CEA, E´cole Polytechnique, F-91128 Palaiseau, France 共Received 8 January 1999; revised manuscript received 22 July 1999兲 An exchange-correlation energy functional beyond the local-density approximation 共LDA兲, based on the exchange-correlation kernel of the homogeneous electron gas and originally introduced by Kohn and Sham, is considered for electronic-structure calculations of semiconductors and atoms. Calculations are carried out for diamond, silicon, silicon carbide, and gallium arsenide. The lattice constants and gaps show a small improvement with respect to the LDA results. However, the corresponding corrections to the total energy of the isolated atoms are not large enough to yield a substantial improvement for the cohesive energy of solids, which remains hence overestimated as in the LDA. 关S0163-1829共99兲01940-2兴

I. INTRODUCTION

Calculations based on the Kohn-Sham1 formulation of density-functional theory 共DFT兲 共Ref. 2兲 have become a prominent tool in condensed-matter physics. Current work is dominated by local-density approximation 共LDA兲 studies, in which the exchange-correlation potential is a local function of the density.1,3 Despite its enormous success, the local-density approximation has some shortcomings which motivated the increasing interest in approaches beyond it, such as the generalized gradient approximation 共GGA兲,4,5 the average-density approximation 共ADA兲, and the weighted density approximation 共WDA兲.6 However, at present none of these methods has replaced the simpler local-density scheme, because they do not always yield systematic and consistent improvements with respect to the LDA, and because, except for the GGA, they could not be implemented in a computationally tractable fashion. It is therefore worthwhile to carry out electronic-structure calculations using other approximations for the exchangecorrelation energy functional. Kohn and Sham1 proposed, together with the LDA, a correction to it which is exact up to the second order in the density fluctuations with respect to the average electron density n; hence it is appropriate for weakly inhomogeneous systems. It involves, as a basic ingredient, the exchange correlation kernel of the homogeheg ជ ជ (r ⫺r ⬘ ), describing the change of the neous electron gas, K xc exchange-correlation potential at rជ induced by a density change at rជ ⬘ . This functional was ruled out since Gunnarson et al.6 showed that different choices of how to include higher-order terms may lead to very different results for the total energy in strongly inhomogeneous systems. In fact, this functional has never been applied to real solids, while the LDA became the state-of-the-art of totalenergy and band-structure calculations. Keeping into account heg based on quanthat several realistic parametrizations of K xc tum Monte Carlo 共QMC兲 calculations are now available, we 0163-1829/99/60共16兲/11329共7兲/$15.00

PRB 60

think that the relative simplicity of this non-LDA 共NLDA兲 functional, and its avoiding gradient expansions to treat inhomogeneity, suggest a further attempt to apply it to solids. In Sec. II we describe in detail the NLDA functional and we propose a derivation which eliminates some arbitrariness in the treatment of higher-order density fluctuations. The choice which we adopt avoids divergencies in the case of strongly inhomogeneous systems. In Sec. III we describe the computational details. In Sec. IV we discuss the results obtained for bulk silicon, diamond, silicon carbide, and gallium arsenide, as well as for atoms and pseudoatoms, and for a surface. The conclusions are drawn in Sec. V. II. NLDA EXCHANGE-CORRELATION ENERGY FUNCTIONAL

In this section we derive a correction to the LDA exchange-correlation 共XC兲 energy for a weakly inhomogeneous system. Although the final result will be coincident with a NLDA XC functional introduced by Kohn and Sham,1 the derivation given here is different. It suggests a wider validity range of the resulting functional, and renders the choice of how to include higher orders less arbitrary. Let us consider the exchange-correlation potential V xc(rជ ) in a weakly inhomogeneous system. We can write LDA „n 共 rជ 兲 …⫹ V xc共关 n 兴 ,rជ 兲 ⫽V xc



heg ជ ជ ˜ ជ ជ drជ ⬘ K xc „r ⫺r ⬘ ;n 共 r ,r ⬘ 兲 …

⫻ 关 n 共 rជ ⬘ 兲 ⫺n 共 rជ 兲兴 ⫹•••,

共1兲

heg ជ ជ ˜ ជ ជ where K xc „r ⫺r ⬘ ;n (r ,r ⬘ )… is the XC kernel 共that is, the functional derivative of the XC potential at rជ with respect to the density at rជ ⬘ ) for the jellium 关homogeneous electron gas 共HEG兲兴 of density ˜n (rជ ,rជ ⬘ ). The first term on the right-hand side of Eq. 共1兲 is obtained assuming, as in the LDA, that the electron density n(rជ ⬘ ) is constant, equal to n(rជ ). The second term is an expansion of V xc( 关 n 兴 ,rជ ) to first order in 关 n(rជ ⬘ )

11 329

©1999 The American Physical Society

11 330

MAURIZIA PALUMMO et al.

⫺n(rជ)兴, while higher-order terms are neglected. Strictly speaking, ˜n (rជ ,rជ ⬘ ) should be taken equal to n(rជ ). However, any choice in the range between n(rជ ) and n(rជ ⬘ ) will not alter the resulting XC potential to first order in n(rជ )⫺n(rជ ⬘ ). The choice of ˜n (rជ ,rជ ⬘ ) will be discussed below. The next problem to solve is to find the XC-energy functional E xc关 n 兴 whose functional derivative ␦ E xc / ␦ n(r) yields the XC potential 共1兲. If ˜n (rជ ,rជ ⬘ ) is a symmetric function of rជ and rជ ⬘ , we can easily show that the functional LDA E xc关 n 兴 ⫽E xc 关n兴⫺





1 4



drជ

heg ជ ជ ˜ ជ ជ drជ ⬘ K xc „r ⫺r ⬘ ;n 共 r ,r ⬘ 兲 …关 n 共 rជ ⬘ 兲 ⫺n 共 rជ 兲兴 2

共2兲 yields the following XC potential: LDA „n 共 rជ 兲 …⫹ V xc共关 n 兴 ,rជ 兲 ⫽V xc



heg ជ ជ ˜ ជ ជ drជ ⬘ K xc „r ⫺r ⬘ ;n 共 r ,r ⬘ 兲 …

⫻ 关 n 共 rជ ⬘ 兲 ⫺n 共 rជ 兲兴 ⫺ 共 1/4兲



drជ ⬙







␦˜n 共 rជ ⬙ ,rជ ⬘ 兲 关 n 共 rជ ⬘ 兲 ⫺n 共 rជ ⬙ 兲兴 2 . ␦ n 共 rជ 兲

drជ ⬘

heg ជ ˜ 共 rជ ⬙ ,rជ ⬘ 兲 … dK xc „r ⬙ ⫺rជ ⬘ ;n

dn 共3兲

The XC-correlation potential 共3兲 is coincident with our starting point 共1兲 up to terms of first order in the density variation n(rជ ⬘ )⫺n(rជ ). Hence the XC-energy functional 共2兲 is consistent up to the second order in 关 n(rជ ⬘ )⫺n(rជ ) 兴 with the XC potential 共1兲. The exchange-correlation energy functional 共2兲 becomes equal to the NLDA of Kohn and Sham if the density arguheg ˜ ជ ជ , n (r ,r ⬘ ), is replaced by the average electron ment in K xc density n. Previous derivations of it were based on a partial summation of the gradient expansion,1 or on assuming weak density fluctuations around the average density.7 The present derivation is based on less restrictive assumptions: considering in fact that the XC kernel is generally of short range in 兩 rជ ⫺rជ ⬘ 兩 共according to Hubbard’s expression,8 we expect that it decays after a few reciprocal Fermi wave vectors兲, the density variation 关 n(rជ ⬘ )⫺n(rជ ) 兴 must be small within this range, in order to allow the neglect of higher-order terms. We refer to this situation as to weak inhomogeneity on the heg ជ ជ (r ⫺r ⬘ ;n), or shortly as scale of the nonlocality range of K xc to weak inhomogeneity. If we consider an inhomogeneous system, the density argument in the HEG XC kernel in Eq. 共2兲 has some arbitrariness: in fact, since the expression is valid up to the second order in density fluctuations, the XC kernel has to be correct to zeroth order in it. This means that in principle any density close to n(rជ ) and n(rជ ⬘ ) can be used therein. The simplest choices which preserve the symmetry in rជ and rជ ⬘ are 共i兲

PRB 60

n 关 (rជ ⫹rជ ⬘ )/2兴 and 共ii兲 关 n(rជ )⫹n(rជ ⬘ ) 兴 /2. Choices 共i兲 and 共ii兲 have been discussed in Ref. 6, where it has been demonstrated that they can lead to very different results in the cases of strongly inhomogeneous systems, where choice 共i兲 can even lead to infinite values of the XC energy. Choice 共i兲 is in fact ruled out by our derivation, since n 关 (rជ ⫹rជ ⬘ )/2兴 has no reason to lie in the density range of n(rជ ) and n(rជ ⬘ ). We show below that choice 共ii兲 is not only consistent with our derivation, but also minimizes the error in real systems with respect to choice 共i兲. heg ជ ជ (r ⫺r ⬘ ;n) becomes The main source of error is that K xc of long range when n is vanishingly small. This is understood on the grounds of self-interaction arguments6 and ocheg ជ ជ (r ⫺r ⬘ ;n), namely curs even in the simplest model of K xc 8 Hubbard’s model. No relevant error occurs in Eq. 共2兲 when both rជ and rជ ⬘ are in regions of small electron density, since in this case the factor 关 n(rជ ⬘ )⫺n(rជ ) 兴 2 vanishes. Serious errors can instead occur when rជ 共or rជ ⬘ ) is in a region of small density, while rជ ⬘ 共or rជ ) is in a region of relevant density. If we make choice 共i兲, n 关 (rជ ⫹rជ ⬘ )/2兴 may happen to be also small, and the resulting kernel of long range. A very big wrong contribution may be added to the XC energy in this case. If, on the other hand, we make choice 共ii兲, the density appearing in the kernel, 关 n(rជ )⫹n(rជ ⬘ ) 兴 /2, is nonvanishing heg ជ ជ (r ⫺r ⬘ ;n) is of short range in all relevant cases. and K xc We adopt hence choice 共ii兲 in the functional 共2兲. The residual error can be checked a posteriori by comparing the expectation values over the electron wave functions of the XC potentials of Eqs. 共1兲 and 共3兲 „with 关 n(rជ )⫹n(rជ ⬘ ) 兴 /2 in place of ˜n (rជ ,rជ ⬘ )…: if the assumption of weak inhomogeneity is valid, they must be very close to each other. We verified that this is not the case for atoms, but it is true for the analyzed semiconductors. The underlying reason for this is that electrons in valence or lower conduction states, which are of interest in most electronic-structure calculations, have very small probabilities of being found in regions of small charge density, where wrong contributions to the XC functional can be generated as described above. In any case we used form 共3兲, which is fully consistent with our ansatz 共2兲 for the XC energy. heg (q;n) of the HEG 共Fourier transThe XC kernel K xc formed from rជ ⫺rជ ⬘ to qជ ) is strictly related to the so-called heg (q;n)⫽ static local-field factor G(q;n): K xc 2 2 ⫺(4 ␲ e /q )G(q;n). We have carried out numerical calculations using different forms of G(q;n): 共i兲 a Hubbard-like 2 (n) 兴 , where q TF(n) is the form, i.e., G(q;n)⫽q 2 /2关 q 2 ⫹q TF Thomas-Fermi wave vector at density n; 共ii兲 two more realistic models for G(q;n), the first based on the known asymptotic behaviors of G(q;n) for small and large q, and on approximate calculations in between—the model of Ichimaru and Utsumi9 共IU兲—and the second based on a parametrization of the QMC data of Moroni and Senatore,10 recently proposed by some of us 共CPOD兲.11 The last two models of G(q;n) differ mainly in the asymptotic behavior for large q 关which in the IU model fails to reproduce the exact result, i.e., G(q;n)⬃q 2 兴, and for the presence, in the IU model, of a logarithmic singularity for

PRB 60

NONLOCAL DENSITY SCHEME FOR ELECTRONIC- . . .

FIG. 1. XC kernels of the homogeneous electron gas, K xc(R), according to the model of Ichimaru and Utsumi 共Ref. 9, full line兲 and to the model of Ref. 11 based on QMC data of Ref. 10 共dashed line兲, as a function of R⫽rk F .

q⫽2q f (q f is the Fermi wave vector兲, which does not appear in the CPOD G(q;n). 共For a more detailed discussion of these differences, see Ref. 11.兲 The resulting XC kernels in real space are quite similar, as shown in Fig. 1, where they are plotted for r s ⫽2. In the first case, the presence of oscillations due to the logarithmic singularity leads to a more difficult convergence of the integrals, without affecting very much the final results. Hence, in the following we adopt the CPOD parametrization; results based on the other two models will be mentioned occasionally for the sake of comparison. III. COMPUTATIONAL DETAILS

Self-consistent electronic-structure calculations were carried out using pseudopotentials and a plane-wave basis. We generated the pseudopotentials on a numerical grid, using the method proposed by Troullier and Martins12 for silicon and diamond, while for Ga and As we followed the scheme introduced by Hamann.13,14 The pseudopotentials have been generated by carrying out atomic calculations within the same NLDA scheme as that used in the solid calculation. All the pseudopotentials were used in the fully separable representation of Kleinman and Bylander,15 after having checked that no ghost state is generated. In the case of GaAs, nonlinear core corrections were included.16 Plane-wave basis sets with a cutoff of 18 Ry for silicon, 20 Ry for GaAs, and 40 Ry for diamond and silicon carbide are needed to achieve good convergence. Ten special points17 were used to sample the irreducible Brillouin zone. The NLDA XC functional 共2兲 and the XC potentials 共1兲

11 331

FIG. 2. Convergence of NLDA correction as a function of r c in bulk silicon, using the model of Ref. 9 共full line兲 and Ref. 11 共dashed line兲.

and 共3兲 were calculated by integration in real space 共on a grid of 16⫻16⫻16 points for Si, GaAs, and diamond, and of 24⫻24⫻24 points for SiC兲, with a cutoff r c in 兩 rជ ⫺rជ ⬘ 兩 . Figure 2 shows the convergence of the total energy of Si 共we plot the correction to the LDA value兲, as a function of r c , for heg . The results are compared with the CPOD form of K xc heg . They those obtained using the UI parametrization of K xc are very similar, but the integrals calculated with the CPOD model converge faster with respect to r c . Clearly the computational effort for the integration increases linearly with the number of grid points used to calculate the integrals. However, two factors contribute in limiting the numerical effort requested by the present NLDA method: 共i兲 the convergence with r c is fast, showing the short-range character of heg both for silicon and for the other materials considered; K xc 共ii兲 since, as we verified, self-consistency effects are very small, all the calculations can be done starting from the LDA wave functions and converge to the true NLDA ground state in a very small number of iterations. We have checked the validity of the assumption of weak inhomogeneity, by calculating, in the case of silicon and diamond, the expectation values of both the XC potentials 共1兲 and 共3兲: they differ by less than 0.01 eV, showing that terms of second order in the density variations are negligible in the XC potential, and therefore that Si and diamond are weakly inhomogeneous systems in the sense discussed in Sec. II. The same is true for cubic SiC, since the energy changes between NLDA and LDA are of the same order as in Si. However in the case of GaAs, where the effect of the large core charge density is partially taken into account through the nonlinear core corrections, the differences are as large as

MAURIZIA PALUMMO et al.

11 332

TABLE I. Ground-state properties of bulk silicon, diamond, silicon carbide, and gallium arsenide as obtained in LDA and NLDA calculations. The GGA data taken from Ref. 5 are reported for completeness. In all cases the pseudopotentials have been generated within the same scheme used for the solid calculation. a 0 (a.u.)

B 0 (Mbar)

E b (eV)

Si

LDA NLDA GGA Expt.

10.17 10.19 10.32 10.26

0.98 0.95 0.87 0.99

5.31 5.28 4.64 4.63

C

LDA NLDA GGA Expt.

6.73 6.75 6.76 6.74

4.51 3.93 4.08 4.42

8.63 8.46 7.80 7.37

SiC

LDA NLDA Expt.

8.15 8.16 8.24

2.25 2.00 2.3

7.42 7.35 6.34

GaAs

LDA NLDA GGA Expt.

10.55 10.60 10.79 10.68

0.77 0.85 0.64 0.75

4.00 4.01 3.31 3.26

0.3 eV, which suggests that the density argument should be chosen carefully. IV. RESULTS A. Bulk systems

In Table I we show the lattice constants (a 0 ), bulk moduli (B 0 ), and binding energies (E b ) obtained for the materials examined: it is evident that all relevant quantities change very little with respect to the LDA. The equilibrium lattice constants increase slightly, while the bulk modulus decreases for Si, C, and SiC and increases for GaAs. Despite the fact that E b remains overestimated as in the LDA, the small changes of the ground-state properties obtained are somehow comforting. In fact, gradient corrections sometimes overcorrect the LDA lattice constants,18,19 as shown also in a recent paper by Fuchs et al.,5 where a new set of GGA results for different solid compounds is reported, and the role of corevalence exchange correlation is investigated. In the present

PRB 60

NLDA scheme, the ground-state properties of the solids remain very close to the LDA results, avoiding overcorrections of the lattice constants which can occur in the case of GGA calculations when the pseudopotentials are consistently generated within the XC scheme used.5 The first column of Table II shows the NLDA and GGA ps ). corrections to the ground-state total energy of solids ( ␦ E sol ps The NLDA E sol is always decreased with respect to the LDA. For silicon, this energy lowering agrees with QMC results obtained by Fahy et al. in Ref. 20, while for diamond, QMC yields a total energy slightly less negative than the LDA one, at variance with our results.20 Anyway, in view of the uncertainty still present in the total energy obtained from different QMC calculations 共mostly due to the usage of a pseudopotential generated within the DFT-LDA and to the variational character of the wave functions; see, for example, the results shown in Refs. 20–22兲, we believe that the discrepancies of the order of a few tenths of an eV per atom are acceptable. For all the materials considered, the NLDA electronic structures 共calculated using NLDA pseudopotentials but at fixed LDA lattice parameter in order to disentangle the indirect effect of the latter兲 show slight differences with respect to the LDA band structure, with a generally very small increase of the gaps. The values of the main gaps with both the LDA and NLDA schemes are shown in Table III. For the silicon case we also calculated the NLDA total heg ជ ជ (r ⫺r ⬘ ,n). It is clear energy using Hubbard’s model of K xc from Eq. 共2兲 that, in this case, the XC NLDA functional is heg is always negative. positively definite, since Hubbard’s K xc In fact, we found an increase of the total energy per atom of about 1.5 eV. On the other hand, it is reasonable that, if the total energy increases, the gap between the filled and the empty states closes. This is indeed our result: the direct gap using Hubbard’s K xc is 0.6 eV smaller than the LDA value. These very different results show that it is very important to use a good description of K xc . B. Atoms

In order to obtain the cohesive energy of solids, one also has to calculate the atomic energies using the same XC functional. This is a stringent test, since the condition of weak inhomogeneity is hardly fulfilled in atoms. We start carrying out all-electron calculations for a number of atoms. In Table IV we compare the results obtained for the all-electron

TABLE II. List of total energy and cohesive energy changes with respect to the LDA results, due to the use of NLDA scheme for the considered materials 共the LDA spin correction for the total energy of pseudoatoms is included兲. ps ␦ E sol (eV)

␦ E atps (eV)

␦ E b (eV)

Si 共NLDA兲 Si 共GGA from Ref. 5兲

⫺0.28 0.60

⫺0.31 0.03

⫺0.03 ⫺0.57

C 共NLDA兲 C 共GGA from Ref. 5兲

⫺0.54 0.31

⫺0.70 ⫺0.71

⫺0.16 ⫺1.02

SiC 共NLDA兲

⫺0.59

⫺0.31 共Si兲; ⫺0.70 共C兲

⫺0.07

GaAs 共NLDA兲

⫺1.03

⫺1.72 共Ga兲; ⫺0.36 共As兲

⫹0.01

NONLOCAL DENSITY SCHEME FOR ELECTRONIC- . . .

PRB 60

11 333

TABLE III. Main direct gaps 共eV兲 for bulk silicon, diamond, and gallium arsenide as obtained in LDA and NLDA calculations at fixed lattice parameter (a 0 ⫽10.16, 6.74, and 10.6 a.u., respectively兲. LDA

NLDA

Si

⌫ X L

2.57 3.57 2.86

2.60 3.63 2.90

C

⌫ X L

5.62 11.21 11.32

5.70 11.32 11.39

GaAs

⌫ X L

0.39 3.99 2.05

0.38 4.01 2.06

ground-state energy of some neutral atoms in the LDA and NLDA schemes with the ‘‘experimental’’ values based on the compilation of Veillard and Clementi 共Ref. 23 as reported in Ref. 24兲, and with recent GGA results.25 The NLDA correction increases with the size of the atoms, but corresponds only to a fraction of the GGA result 共which reproduces well the experimental data兲. The reason for this partial failure is evident in Fig. 3, where the NLDA V xc potential of the neon atom is shown together with the LDA and GGA ones. The differences are clear: there is a region, at intermediate distances from the nucleus, where the NLDA XC potential becomes more negative than the LDA potential, but at a large distance it remains substantially equal, with no correction to the wrong 共exponential兲 decay of the LDA V xc(r). The GGA has this shortcoming too, but compensates for it in the region near the GGA is more attractive than ours and the nucleus, where V xc GGA LDA V xc . Whether or not V xc is closer to the true exchange-correlation potential than ours in this region is not clear at present. It is possible that the good agreement between GGA and ‘‘exact’’ total energies of atoms is a consequence of the cancellation of two errors, namely the lack of the Coulombic tail at long range and a possible overactive behavior at small distances, the reciprocal cancellation of TABLE IV. All-electron ground-state total energy 共Hartree兲 of some atoms from He to argon, obtained in the LDA,GGA 共Ref. 25兲 and NLDA schemes, and compared with the experimental values 共from Ref. 23, as reported in Ref. 24兲. Each calculated total energy includes the spin polarization corrections obtained in the LDA scheme. g.s. energy 共a.u.兲 4

Be C 8 O 10 Ne 11 Na 14 Si 18 Ar 6

LDA

NLDA

GGA

Expt.

⫺14.44 ⫺37.46 ⫺74.50 ⫺128.18 ⫺161.38 ⫺288.09 ⫺525.65

⫺14.51 ⫺37.62 ⫺74.81 ⫺128.77 ⫺162.11 ⫺289.23 ⫺527.59

⫺14.64 ⫺37.78 ⫺74.99 ⫺128.94 ⫺162.25 ⫺289.33 ⫺527.54

⫺14.67 ⫺37.84 ⫺75.07 ⫺128.94 ⫺162.25 ⫺289.34 ⫺527.54

FIG. 3. The XC LDA potential 共full line兲 compared with the potential obtained including self-consistently the NLDA correction 共dashed line兲 and with the GGA potential 共dotted line兲 for the Ne atom.

these two errors being optimized by the suitable choice of parameters made within the GGAs. For the calculation of cohesive energies, however, the allelectron results for the atoms must be completed with the calculation of the atomic energies within the same pseudopotential approach used for the solid. Hence, we generated the pseudopotentials including the NLDA exchangecorrelation energy and applied the NLDA functional to pseudoatoms, in order to have a consistent approach in the atomic and solid calculations. In Fig. 4 we report the silicon pseudopotentials, within the LDA and the NLDA scheme: the curves obtained in the two schemes look quite similar, the main differences being confined to the core region (r ⬍1 a.u.). In this region, which is not really important for the computation of matrix elements, being weighted by the r 2 dr volume element, the NLDA pseudopotentials are slightly less attractive than the LDA ones. Less evident, but more important, are the differences for r⬎1 a.u., where the NLDA curves are slightly deeper than the LDA ones. ps between the ground-state total enerThe differences ␦ E at gies obtained in the NLDA and LDA schemes for the pseudoatoms considered here are reported in the second column of Table II. The corresponding differences, taken from Ref. 5, between the GGA and the LDA schemes are also reported for comparison. In our NLDA calculations, the atoms are always treated as spin-unpolarized systems, assuming that any contribution of spin polarization is well described at the local-spin density 共LSD兲 level. 共This conjecture is substantiated by Fuchs et al.’s calculations,5 where spin-polarization energies calculated within GGA and LSD differ by about 0.1 eV.兲 Contrary to the all-electron

MAURIZIA PALUMMO et al.

11 334

PRB 60

DFT-LDA calculation, the buckled dimer 2⫻1 structure is found to be lower in energy than the symmetric dimer one by 0.2 eV/dimer, while the p(2⫻2) and c(4⫻2) reconstructions 共which differ only by 3 meV per dimer兲 show a decrease of total energy with respect to the buckled 2⫻1 phase of 0.02 eV/dimer.28 We have done a slab calculation with twelve atomic layers, fouro vacuum layers, and four special points to sample the irreducible part of the Brillouin zone, using both the LDA and the NLDA scheme. We obtained that the total energy of the slab, in the NLDA scheme, is 0.12 eV/atom lower with respect to the LDA one, which is a decrease of the same order as the one obtained in bulk silicon. However, the equilibrium geometry does not show any appreciable difference with respect to that obtained using the LDA exchange and correlation potential. V. CONCLUSIONS

FIG. 4. l⫽0, l⫽1, and l⫽2 components of the unscreened pseudopotentials for silicon within the NLDA 共dashed lines兲 and LDA 共full lines兲.

case, little improvement with respect to the LDA is obtained for the total energies, since most of the NLDA effect shown in Table IV originates within the core region. Hence the calculated differences in the binding energy ␦ E b are small, and the values of the NLDA E b , contrary to the case of the GGA, do not improve substantially with respect to the LDA values 共see also the last column of Table I兲. Those calculated within the GGA, indeed, correct most of the error of the LDA. C. Surfaces

Although the inclusion of the NLDA functional introduced here does not show strong differences with respect to the LDA for the two limit cases, bulk and atomic systems, for the sake of completeness we also report briefly the results obtained for a surface. We calculated the ground-state energy and the equilibrium structure of the (2⫻1) reconstruction of Si共100兲 using both the standard LDA and the NLDA functional. This surface is extensively studied both from the experimental and theoretical points of view, due to its technological importance.26 Now, it is generally accepted that its atomic structure is characterized by the presence of asym¯ 兴 direction, metric buckled surface dimers along the 关 011 which implies at least a (2⫻1) reconstruction. Two other energetically competing 共favored兲 reconstructions are possible, within the same top-atom bonding characteristic. The alternation of dimer buckling along the 关 011兴 dimer rows leads to a p(2⫻2) reconstruction, while the alternation of buckling angles also along the direction perpendicular to the rows leads to the c(4⫻2) phase.27 In a well-converged

An old non-LDA exchange-correlation functional, originally derived by Kohn and Sham for weakly inhomogeneous systems, has been implemented for usage in electronicstructure calculations. A new derivation, and a well defined treatment of higher-order fluctuations, have been devised, suggesting that this functional should be reliable for systems with slow, rather than weak, density variations 共i.e., with density variations smooth on the scale of the inverse Fermi wave vector兲. NLDA results for the lattice constant of Si, diamond, cubic SiC, and GaAs are slightly better than the LDA ones when compared with experiment. The overcorrection which sometimes occur in the case of the generalized gradient approximations is avoided. Negligible changes with respect to LDA have been found also for a surface, Si(100)2⫻1. On the other hand, the application of this XC functional to atoms 共pseudoatoms兲 improves the total energies only partially 共slightly兲 with respect to the LDA values. Hence, no substantial improvement comes out for the cohesive energy of solids. The reason for this failure is probably due to the fact that the Coulombic long-range tail of V xc , which is important in the case of atoms, is still lacking. In general, we obtain changes with respect to LDA which are smaller than those needed to match the experimental values. This is a consequence of the conservative ansatz made about the density argument in Eqs. 共1兲 and 共2兲, meant to avoid long-range components of K xc and possible divergences of the XC energy. Even though this goal has been obtained, the NLDA part of the XC potential is generally underestimated. This leads to important discrepancies with experiments in the case of atoms, whose binding energies are not fully corrected with respect to the LDA. The GGA, instead, yields an almost complete correction. However, in some cases, e.g., for bulk semiconductors, the present description may be closer to reality. In conclusion, we believe it is worthwhile to continue the search for density functionals of the type considered here, by trying to give a better description of the long-range tail of the exchange-correlation kernel K xc . ACKNOWLEDGMENTS

We acknowledge useful discussions with R. Resta and G. Senatore in the early stage of this work. This work was sup-

PRB 60

NONLOCAL DENSITY SCHEME FOR ELECTRONIC- . . .

11 335

ported by the INFM Parallel Computing Initiative, which provided us with computer resources on an SGI Origin 2000 system at CINECA 共Interuniversity Consortium of Northeastern Italy for Automatic Computing兲. We are grateful to

S. Goedecker for his efficient code for fast Fourier transforms.29,30 This work was supported in part by the European Community program ‘‘Human Capital and Mobility’’ 共Contract No. ERB CHRX CT930337兲.

W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 共1965兲. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 共1964兲. 3 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 共1980兲; J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 共1981兲. 4 J. P. Perdew, Phys. Rev. Lett. 55, 1665 共1985兲; A. D. Becke, Phys. Rev. A 38, 3098 共1988兲; D. C. Langreth and M. J. Mehl, Phys. Rev. Lett. 47, 446 共1981兲; J. P. Perdew, K. Burke, and M. Ernzerhof, ibid. 77, 3865 共1997兲. 5 M. Fuchs, M. Bockstedte, E. Pehlke, and M. Scheffler, Phys. Rev. B 57, 2134 共1998兲, and references therein. 6 O. Gunnarson, M. Johnson, and B. I. Lundqvist, Phys. Rev. B 20, 3136 共1979兲. 7 L. J. Sham, Phys. Rev. B 7, 4357 共1973兲. 8 J. Hubbard, Proc. R. Soc. London, Ser. A 243, 336 共1958兲. 9 S. Ichimaru and K. Utsumi, Phys. Rev. B 24, 7385 共1981兲. 10 S. Moroni, D. M. Ceperley, and G. Senatore, Phys. Rev. Lett. 75, 689 共1995兲. 11 M. Corradini, R. Del Sole, G. Onida, and M. Palummo, Phys. Rev. B 57, 14 569 共1998兲. 12 N. Trouiller and J. L. Martins, Phys. Rev. B 43, 1993 共1991兲. 13 D. R. Hamann, Phys. Rev. B 40, 2980 共1989兲. 14 M. Fuchs and M. Scheffler, Comput. Phys. Commun. 119, 67 共1999兲. 15 L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 共1982兲.

16

1 2

S. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 共1982兲. 17 H. J. Monkorst and J. D. Pack, Phys. Rev. B 13, 5188 共1976兲. 18 A. Garcia, C. Elsasser, J. Zhu, and S. G. Louie, Phys. Rev. B 46, 9829 共1992兲, and references therein. 19 A. Dal Corso, S. Baroni, and R. Resta, Phys. Rev. B 49, 5323 共1984兲. 20 S. Fahy, X. W. Wang, and S. G. Louie, Phys. Rev. B 42, 3503 共1990兲. 21 X. P. Li, D. M. Ceperley and R. M. Martin, Phys. Rev. B 44, 10 929 共1991兲. 22 B. Kralik, P. Delaney, and S. Louie, Phys. Rev. Lett. 80, 4253 共1998兲. 23 A. Veillard and E. Clementi, J. Chem. Phys. 49, 2415 共1968兲. 24 M. T. Carroll, R. F. W. Bader, and S. H. Vosko, J. Phys. B 20, 3599 共1987兲. 25 I. Lee and R. Martin, Phys. Rev. B 56, 7197 共1997兲. 26 M. Palummo, G. Onida, R. Del Sole, and B. Mendoza, Phys. Rev. B 60, 2522 共1999兲. 27 J. Northrup, Phys. Rev. B 47, 10 032 共1993兲. 28 A. Ramstad, G. Brooks, and P. J. Kelly, Phys. Rev. B 51, 14 504 共1995兲. 29 S. Goedecker, Comput. Phys. Commun. 76, 294 共1993兲. 30 S. Goedecker, SIAM 共Soc. Ind. Appl. Math.兲 J. Sci. Stat. Comput. 18, 1605 共1997兲.