text:
# Comparative study of density functional theories of the exchange-correlation hole and energy in silicon

## I Introduction

Density functional theory (DFT) is the leading theoretical tool for the study of the material properties of solids. It is based on the characterization of the ground-state energy of the inhomogeneous electron gas as a functional of the density, whose optimization condition can be expressed in terms of the self-consistent single-particle Kohn-Sham equation.  The success of the method lies in the ease of use and surprising effectiveness of the standard approximation, the local density approximation (LDA). This models the key component of the functional, the exchange-correlation energy $`E_{xc}[n]`$, and the exchange-correlation potential used in the Kohn-Sham equation with the assumption that the inhomogeneous electron gas at any location in space behaves like the homogeneous electron gas at a density equal to the local density. As a result, the formidably complex and nonlocal relationship between the energy and the density due to interparticle correlations in the inhomogeneous environment can be converted into a local functional dependent on input from the comparatively simpler homogeneous electron gas. This approximation has been quite successful at obtaining useful and suprisingly accurate estimates for bulk ground-state properties of many solids and the basic ground-state structure of solids, surfaces and molecules.

However, there remains a significant need for the development of more accurate density functionals in several areas. Applications in molecular and solid-state chemistry require the calculation of total energies to at least an order of magnitude greater precision than the LDA can provide. Excited state properties, in particular the band-gap of semiconductors, are not typically obtainable with quantitative accuracy. In addition, in materials in which electron correlations are important, the LDA often fails to predict qualitative ground-state properties. There have been numerous attempts in recent years to develop exchange-correlation functionals beyond the LDA, including generalized gradient approximations (GGA’s)  based on gradient expansions about the LDA, model nonlocal potentials,  self-interaction corrections,  hybrids with Hartree-Fock  and configuration interaction  or other many-body theories,  and orbital-dependent functionals.  These have led to significant improvements in accuracy of DFT calculations; but a functional of systematically quantitative accuracy for many quantum chemistry and solid-state applications has yet to be achieved.

A key relation that has been particularly fruitful in the development of density functional theory is the adiabatic connection formula  which relates the exchange-correlation energy to a coupling-constant integrated exchange-correlation hole. The exchange-correlation hole, $`n_{xc}`$, measures the change in density from its mean value at each point in a system, given the observation of an electron at one particular point, providing a simple visualization of the effects of electron correlation in an inhomogeneous material. By integrating this quantity over a family of systems characterized by interaction coupling constant $`\lambda `$ and kept at fixed density, the kinetic energy cost to create the hole at full coupling, $`\lambda =1`$, is incorporated into it. As a result, the Coulomb interaction energy of the integrated hole reflects the total energy associated with creating the hole.

The adiabatic connection formula, by making explicit the relation between the exchange-correlation energy and the exchange-correlation hole has been the impetus, both explicitly and implicitly, to many proposed corrections to the LDA. The average density approximation  (ADA) and weighted density approximation  (WDA) make explicit use of this relation to construct nonperturbative nonlocal density functionals. WDA in particular has attracted attention as a potential tool with applications for quantum chemistry,  metallic surfaces,  and bulk solid-state structure.  It is apparently limited by a lack of consistency in its behavior  and typically requires the use of various context-dependent models or extensions to achieve optimal results for different applications.  The PW91  and PBE  forms of the GGA explicitly build in a number of exact and approximate properties of the system average of the exchange-correlation hole and energy using a gradient expansion of $`n_{xc}`$  as a framework. Hybrid models,  which have been particularly successful in achieving quantitative improvements of binding and total energies in molecules, may also be motivated as approximations to the coupling constant integration of the adiabatic connection formula.

In developing many of these approaches, many-body calculations of the exchange-correlation hole and related quantities in real systems, particularly in atoms and simple molecules, have played a useful role, helping to confirm assumptions of approximate functionals or point out specific areas for improvement.  There have been few such studies for real solids, however.

A second fruitful decomposition of $`E_{xc}`$ has been to separate it into exchange and correlation components. The exchange energy, which contains the correlation due to Fermi statistics that occurs for the noninteracting or $`\lambda =0`$ limit, may be explicitly written in terms of the occupied single-particle orbitals obtained in this limit, and thus it and its functional derivative with density can in principle be obtained exactly and self consistently. Applications along this line include the optimized effective potential (OEP) method,  an implementation of the exact exchange potential which has been limited in practice to free atoms and other simple geometries, and the method of Krieger, Li and Iafrate (KLI)  which provides an approximate evaluation of the exchange functional by means of a simplified single-particle Green’s function. Recently, an exact-exchange formalism applicable to solids has been developed by Städele et al.  These methods are of particular interest because they are self-interaction free and are promising candidates for the quantitative treatment of excited states.

On the other hand, these advances expose an important problem: many traditional methods such as LDA rely for their success on the close relation between exchange and correlation, in which the correlation hole tends to correct for much of the nonlocal behavior in the exchange hole.  As a result, the error for the sum of exchange and correlation energies is typically significantly smaller than for either separately.  It is of interest, therefore, to study how the correlation hole alone behaves in real materials and how current functionals for correlation fare with respect to accurate many-body calculations. Or from the complementary point of view, it becomes necessary to consider to what extent modeling orbital effects in the correlation hole become important once they are included in exchange.

Variational Quantum Monte Carlo (VMC) calculations of the coupling-constant pair correlation function for the Si crystal have been reported previously,  including an analysis in depth of the exchange-correlation energy density obtained by the adiabatic connection theory. In this work we discuss in more detail the exchange-correlation hole, its decomposition into exchange and correlation and the relation between the angle-averaged hole and the energy density. We present analysis along two lines: first a detailed comparison of the VMC exchange-correlation hole and that of several DFT models as a function of position in the crystal, in order to search for systematic trends in how these functionals describe the exchange-correlation hole of a prototypical semiconductor material. Secondly, we look at how the exchange and correlation holes separately behave in the bonding region of Si, and observe how orbital effects in the exchange generate related effects in the correlation hole.

The paper is organized as follows: Sec. II provides theoretical background on the exchange-correlation hole, and the various models used to describe it. We briefly discuss the method used to obtain computational results in Sec. III. Results for spherically averaged holes are presented in Sec. IV, and a detailed analysis of the correlation and exchange holes in Sec. V. In Sec. VI we discuss issues raised by our data analysis concerning the improved treatment of nonlocal density information in DFT’s, as well as orbital effects in the Si correlation and exchange holes that may be of relevance in improving upon exact exchange methods. Our conclusions are presented in Sec. VII.

## II The exchange-correlation hole in density functional theory

### A Basic theory

We consider a family of systems parameterized by a coupling-constant $`\lambda `$

$$H(\lambda )=\underset{i}{}\left(\frac{_i^2}{2}+V_{ext}(𝐫_𝐢)+V_\lambda (𝐫_𝐢)\right)+\underset{i>j}{}\frac{\lambda }{|𝐫_𝐢𝐫_𝐣|}$$
(1)
where $`V_{ext}`$ is the external potential and $`V_\lambda `$ a potential added to keep the ground-state density invariant. The units here and elsewhere in the paper are in Hartree atomic units unless otherwise indicated. Two limits of interest include the noninteracting system, $`\lambda =0`$, in which $`H`$ reduces to the standard Kohn-Sham equation, with the adjustable potential equal to the Kohn-Sham potential. The second is at $`\lambda =1`$, where $`V_\lambda `$ vanishes and one recovers the original fully interacting many-body system. The ground-state wave function is a Slater determinant of single-particle orbitals in the first case and the true many-body ground state in the second.

The exchange-correlation hole is defined for a given value of $`\lambda `$ as the change in the ground-state expectation of the density at one point in the system $`𝐫+𝐮`$ given the observation of an electron at some other point $`𝐫`$:

$`n_{xc}`$ $`(𝐫,𝐫+𝐮;\lambda )=`$ (3)
$`{\displaystyle \frac{1}{n(𝐫)}}{\displaystyle \underset{i}{}}\delta (𝐫𝐫_i){\displaystyle \underset{ji}{}}\delta (𝐫+𝐮𝐫_j)_\lambda n(𝐫+𝐮).`$
Here $`_\lambda `$ indicates an average over the ground-state wave function of $`H(\lambda )`$ and $`n(𝐫)=_i\delta (𝐫𝐫_i)_\lambda `$ is the $`\lambda `$-invariant ground-state density. A related quantity, the pair correlation function $`g(𝐫,𝐫+𝐮)`$, is a measure of the pair density relative to that expected for uncorrelated electrons with the same density distribution. The exchange-correlation hole is expressed in terms of $`g`$ as

$$n_{xc}(𝐫,𝐫+𝐮;\lambda )=n(𝐫+𝐮)[g(𝐫,𝐫+𝐮,\lambda )1].$$
(5)
The importance of the exchange-correlation hole to density functional theory lies in its connection  to the exchange-correlation energy $`E_{xc}`$:

$$E_{xc}[n]=\frac{1}{2}d^3rn(𝐫)d^3r^{}_0^1𝑑\lambda \frac{n_{xc}(𝐫,𝐫^{},\lambda )}{|𝐫𝐫^{}|}.$$
(6)
$`E_{xc}`$ includes all the contributions to the energy due to correlations, that is, beyond the noninteracting and Hartree energy. The exchange-correlation energy $`E_{xc}[n]`$ is analyzed in density functional theory in terms of the exchange-correlation energy density $`e_{xc}(𝐫)`$ and the exchange-correlation energy per particle $`ϵ_{xc}(𝐫)`$

$$E_{xc}=d^3rn(r)ϵ_{xc}(𝐫)=d^3re_{xc}(𝐫).$$
(7)
Relating this expression to that for $`E_{xc}`$ in terms of the exchange-correlation hole $`n_{xc}`$ one has

$$ϵ_{xc}(𝐫)=\frac{1}{2}d^3u\frac{n_{xc}(𝐫,𝐫+𝐮)}{u},$$
(8)
where $`n_{xc}(𝐫,𝐫+𝐮)`$ is the coupling-constant integrated hole, $`_0^1𝑑\lambda n_{xc}(𝐫,𝐫+𝐮,\lambda )`$. The exchange-correlation energy per particle thus has the natural interpretation as the interaction energy of the particle with its exchange-correlation hole.

The exchange-correlation hole and energy are frequently decomposed according to separate exchange and correlation components. The exchange hole $`n_x`$ is that of the noninteracting or $`\lambda =0`$ ground state \[Eq. (LABEL:eqnxc)\] and contains correlations due solely to Fermi statistics. The remainder, $`n_c`$, of the contributions to the total $`n_{xc}`$ are those induced by introducing pair correlations into the ground-state wave function, representing density fluctuations caused by the Coulomb interaction. The two components satisfy important sum rules which reflect the conservation of particle number:

$`{\displaystyle d^3un_x(𝐫,𝐫+𝐮)}`$ $`=`$ $`1,`$ (9)
$`{\displaystyle d^3un_c(𝐫,𝐫+𝐮)}`$ $`=`$ $`0.`$ (10)
A key to the practical exploitation of the relation between $`n_{xc}`$ and $`E_{xc}`$ is that, since the Coulomb interaction depends only on interparticle distance, much of the information contained in $`n_{xc}`$ is not needed and may be ignored in the development of models. In particular, to obtain the energy per particle $`ϵ_{xc}(𝐫)`$ at a given point $`𝐫`$ the angular variation in $`n_{xc}`$ may be integrated out to leave a function depending only on the distance $`u`$ from $`𝐫`$,

$$n_{xc}(𝐫,u)=\frac{d\mathrm{\Omega }_u}{4\pi }n_{xc}(𝐫,𝐫+𝐮),$$
(11)
with $`ϵ_{xc}(𝐫)`$ in Eq. 8 given in terms of the integral over $`u`$ of the weighted angle-averaged hole, $`2\pi un_{xc}(𝐫,u)`$. The DFT models we study in this paper can be viewed as models of this angle-averaged hole.

### B Models for the exchange-correlation hole

The LDA model for $`n_{xc}`$ may be obtained by replacing the true hole about $`𝐫`$ with that of the homogeneous electron gas at the local density, $`n(𝐫)`$:

$$\text{ }n_{xc}^{LDA}(𝐫,u)\text{ }=n(𝐫)\{g^{heg}[u,r_s(𝐫)]1\}.$$
(12)
The factor $`r_s=(3/4\pi n)^{1/3}`$ is the average distance between electrons in the homogeneous gas of density $`n`$. It provides a useful measure of the correlation length for the exchange-correlation hole, as the width of the exchange hole scales linearly with $`r_s`$ and that of the correlation hole shows similar behavior for the range of densities $`1.4<r_s<4.4`$ valid for the Si crystal valence.

The WDA  attempts to account for the inhomogeneity of the density in realistic systems by incorporating it in the density prefactor of $`n_{xc}`$

$$\text{ }n_{xc}^{WDA}(𝐫,u)\text{ }=\frac{d\mathrm{\Omega }_u}{4\pi }n(𝐫+𝐮)\{g^{model}[u,\overline{r}_s(𝐫)]1\}.$$
(13)
This form of the exchange-correlation hole treats the hole at $`𝐫+𝐮`$ as the modification of the density at that point rather than that at the center of the hole. This is a property of the exact hole that leads to a significant departure from the homogeneous case when the density varies considerably within the effective range of the hole. The major approximation made is the use of a scaling form for the pair correlation function where for each $`𝐫`$ an averaged length scale $`r_s(𝐫)`$ is determined to select $`g`$ from a family of suitably chosen pair correlation functions $`g(u,r_s)`$. This length scale is fixed at each point fixed by the particle sum rule

$$d^3un(𝐫+𝐮)\{g^{model}[u,\overline{r}_s(𝐫)]1\}=1.$$
(14)
A natural choice for $`g(u,r_s)`$ in an extended system is that of the homogeneous electron gas; this form will be considered in this paper. The WDA scaling approximation has the ability to fit certain properties of the true $`n_{xc}`$ unattainable by the LDA; in particular it significantly improves on the description of the exchange-correlation energy per particle of atoms when a single electron is removed, obtaining the correct limiting value of $`1/2R`$ versus distance from the atom $`R`$

The ADA  is another method to incorporate into a model of the exchange-correlation hole information about the nonlocal density variation in the vicinity of the hole. It preserves the homogeneous gas model for $`n_{xc}`$ but determines the $`r_s`$ factor at which it is evaluated from an average of the density in the neighborhood of the reference point:

$$\text{ }n_{xc}^{ADA}(𝐫,u)\text{ }=\overline{n}(𝐫)\{g^{heg}[u,\overline{r}_s(𝐫)]1\}$$
(15)
with

$$\overline{n}(𝐫)=d^3r^{}n(𝐫^{})w(|𝐫𝐫^{}|,\overline{r}_s).$$
(16)
Here $`\overline{r}_s=(3/4\pi \overline{n})^{1/3}`$ and $`w`$ is a normalized weighting factor. The weight $`w`$ is chosen to obtain the correct $`E_{xc}`$ in the linear response limit and has an effective range set by $`\overline{r}_s`$. The averaging procedure can in this sense be understood as the selecting of a physically reasonable correlation length for the hole from an average of the density over the physical range of the hole, particularly of value in cases where the density varies rapidly within this range. Unlike the WDA, the ADA does not try to provide a correction to the shape of the $`n_{xc}`$.

Note that with the use of a radially isotropic $`g`$, only a spherically averaged density enters into the definition of, and sum rule for, the WDA exchange-correlation hole:

$$\text{ }n_{xc}^{WDA}(𝐫,u)\text{ }=n(𝐫,u)\{g^{heg}[u,\overline{r}_s(𝐫)]1\},$$
(17)
and

$$n(𝐫,u)=\frac{d\mathrm{\Omega }_u}{4\pi }n(𝐫+𝐮).$$
(18)
A similar result holds for the ADA, as the density averaging function $`w`$ is a function only of interparticle distance. This averaging reduces the information about the nonlocal density actually used by these methods from the full complexity of the density variation about each reference point to a simpler and often slowly varying radial function $`n(𝐫,u)`$. Knowledge of this function at each point in space, along with the model used for $`g`$ in the WDA and $`w`$ in the ADA, specifies $`E_{xc}`$ within each approximation.

Another line of attack in improving the exchange-correlation hole and energy, the GGA, has been to design corrections to the LDA hole using information about system inhomogeneity contained in the local density gradient. One implementation of this approach  has led to the construction of a GGA model for $`n_{xc}`$ In this and other GGA models the measure of inhomogeneity is the absolute value of the gradient, $`|n(𝐫)|/n(𝐫)`$. Second-order effects proportional to the Laplacian of the density, $`^2n(𝐫)/n(𝐫)`$, are mapped to a gradient form by an integration by parts. This transformation has no effect on system-averaged quantities, namely $`E_{xc}`$ and the average over $`𝐫`$ of $`n_{xc}`$. However, the local quantities $`n_{xc}(𝐫,𝐫+𝐮)`$ and $`ϵ_{xc}(𝐫)`$ in the GGA are as a result no longer defined in terms of the adiabatic formulae of Eqs. \[LABEL:eqnxc\] and  and not directly comparable with our VMC data.

## III System and computational details

A description of the computational method used in obtaining numerical many-body expectations of the pair correlation function and related quantities in Si has been presented in a previous paper.  We present here a short summary of the details relevant to our discussion.

For $`V_{ext}`$ in Eq. (1) we use a norm conserving nonlocal LDA pseudopotential which replaces the electrons in the atom core. The system and expectations involve valence electrons and their correlations. A simulation cell consisting of a $`3\times 3\times 3`$ lattice of fcc primitive cells of the diamond lattice containing 216 valence electrons is used in obtaining expectations. Finite size effects in the electron-electron interaction were taken into account by replacing the Coulomb interaction in Eq. (1) with a truncated potential into which the long-range Hartree interaction is incorporated. This form has proved to be successful in reducing finite size errors in the Coulomb energy as compared to the conventional Ewald interaction.

To describe the ground state for a given value of coupling constant, a Slater-Jastrow type wave function was used, consisting of a product of Slater determinants multiplied by a trial many-body correlation factor:

$$\psi ^\lambda =\mathrm{exp}\left(\underset{ij}{}u^\lambda (|𝐫_i𝐫_j|)+\underset{i}{}\chi ^\lambda (𝐫_i)\right)DD.$$
(19)
The orbitals used in the Slater determinants for up and down spins $`D_\sigma `$ are determined from an LDA calculation and include all the valence-band Bloch orbitals of the primitive cell periodic on the supercell. The Jastrow factor introduces variationally explicit electron-electron correlations with one-body corrections to optimize the density. It is described in detail in the paper by Hood et al.  A variational ground-state wave function optimized over 22 variational parameters was obtained for five different values of the coupling-constant, and recovered 88$`\%`$ of the correlation energy as compared to nearly exact diffusion Monte Carlo calculations for the same system.

The pair correlation function $`g(𝐫,𝐫^{})`$ \[Eq.(5)\] was obtained from the ground-state wave function for each $`\lambda `$ in order to obtain coupling-constant integrated exchange-correlation holes and energies. The pair correlation function was expanded in symmetrized plane-waves, taking full advantage of the symmetries of the pair correlation function in the crystal to reduce the basis set to a reasonable number. Expectations for the coefficients for symmetrized plane-waves were evaluated statistically using the Monte Carlo method. The resulting data are limited by plane-wave cutoff and statistical errors which can be checked by comparing the numerically calculated and exact values of the exchange hole. The resulting variation in $`g_x`$ was between 1% and 6%. The error due to basis set cutoff is most noticeable near the atom core where, despite the pseudopotential treatment, the rapid variation of the density and the valence orbitals leads to large wave-numbers in the exchange hole.

## IV Angle-averaged exchange-correlation hole

Figure 1 shows as a function of interelectron distance $`u`$ the angle-averaged correlation, exchange and exchange-correlation holes, weighted by $`2\pi u`$, for the VMC calculation of Hood et al.  and several DFT models at several representative points in the Si crystal. The pair correlation function of the homogeneous electron gas used to generate the DFT holes is obtained from an accurate model form.  The distance weighting factor $`2\pi u`$ is chosen so that the integral of the plotted function yields the exchange-correlation energy per particle at the reference point. These values are listed separately in Table I.

In order to give an idea of the environment surrounding each reference point and to compare with the WDA and ADA forms for $`n_{xc}`$, the angle-averaged density $`n(𝐫,u)`$, the average density at a distance $`u`$ from the reference point $`𝐫`$, is shown at the top of each plot. The density is to be compared to the average over the unit cell of 0.0296 a.u. The angle-averaged density for each of the four reference points shown indicate four different environments typical of the pseudopotential model Si crystal. In the first, (a), the reference electron is at a density considerably higher than the unit-cell average, and the density decreases significantly over the length-scale of the VMC exchange-correlation hole. The antibond position (b) is typical of intermediate density points where the local density is close to the unit-cell average and the angle-averaged density in the vicinity of the reference point shows only modest variations over long distances. The atom center (c) is typical of the “pathological” situation inside the pseudopotential core. Here the density increases by an order of magnitude, from roughly 10$`\%`$ to 200$`\%`$ of the unit-cell average, over roughly one-third of the Si bond length of 4.44 a.u., before slowly decaying back to the average. Such an environment of dramatic density change over the length scale of the hole indicates the likely failure of the LDA in this region, which has in fact been observed in the comparison of LDA and VMC energy densities.  Finally, the interstitial point (d) shows the generic features at low density outside both valence and core. Here the angle-averaged density also changes considerably from the local density but at a much more gradual rate.

In the high-density case of the bond center, Fig. 1(a), the exchange-correlation hole closely follows that determined by the local density, with $`r_s=1.397`$ a.u. The average density within the effective range of the hole is lower than the local density, leading to a slightly larger $`\overline{r}_s`$ for the ADA and an overall better fit of the hole. The WDA hole \[Eq. (17)\] shows the effect of two competing factors. The value of $`\overline{r}_s`$ used in the scaled pair correlation function $`g^{heg}`$ is significantly larger than the local $`r_s`$. The resulting pair correlation function is wider and deeper than that of the LDA. However, the density prefactor $`n(u)`$ in Eq. (17) drops rather quickly and “quenches” the exchange-correlation hole at larger $`u`$. Eventually, as the density settles to its asymptotic value of the unit cell average, the hole has a long-range tail again dominated by the shape of the pair correlation function. The net result is a considerable overcorrection of the LDA model with respect to the VMC data.

At the antibond point, the hole is wider and less deep than at the bond center, consistent with the variation of the homogeneous electron gas hole with respect to density. Neither form of nonlocal averaging results in a significant change from the LDA form of $`n_{xc}`$ since the angle-averaged density is almost constant in the vicinity of the reference electron. The exchange-correlation energy per particle, determined from the integral of the curves in Fig. 1, has roughly the same modest deviation from the VMC data in each case, as shown in Table I. This leads to a significant difference in the energy density because the density is still large at this point.

The holes at the bond center and antibond point are fairly representative of those at other points of intermediate-to-high density in the Si crystal. The corrections introduced by the nonlocal DFT’s are quite small because of the relatively small variation in the angle-averaged density within the effective range of the hole. These high density corrections to $`n_{xc}`$ and $`e_{xc}`$ have a sensitive dependence on position in the crystal, which appears to lack an obvious relation to the positional trends in the VMC data. As a result, plots of the deviation of the ADA and WDA energy density from the VMC data feature quite complex “dapple” patterns that remain after the worst errors of the LDA model are smoothed out.

For the two low-density cases, however, the nonlocal density averaging techniques have a considerable impact on the shape and quality of the exchange-correlation hole. In these cases the shape of hole is dominated by the rapid increase of the density with $`u`$.

At the atom center, Fig 1(c), the rapid change in density leads to a far more compact exchange-correlation hole than that predicted by the LDA at the low local density ($`r_s`$ = 4.12 a.u.). The minimum in the weighted hole, where it contributes the most to the exchange-correlation energy, is 1.6 a.u. from the reference point, less than half that of the LDA hole. Roughly $`90\%`$ of the total contribution to $`ϵ_{xc}(𝐫)`$ comes from within a radius of 3.0 a.u., significantly less than a bond length (4.44 a.u.). In this case, both ADA and WDA do significantly better at determining the overall depth and length scale of the hole over the energetically important region of 0.5 to 3.5 a.u. In addition, WDA clearly matches the shape of the VMC hole quite closely. As a result, the ADA and WDA obtain fairly accurate estimates of the energy of the exchange-correlation hole as shown in Table I, while LDA obtains only half the VMC value.

The interstitial point, Fig. 1(d), has low-density features qualitatively similar to the atom center but, consistent with a more gradual change in density, is considerably shallower and more spread out. The WDA has a particularly good overall agreement with the VMC hole on a point-by-point basis (although the ADA has a total energy slightly closer to the VMC value). However, possibly due to the more gradual change in density within the vicinity of the hole, the LDA does a much better job of matching the energy of the hole than near the atom center. The large differences in the shape of the hole largely cancel out in the integral so that the energy of the hole is only 20$`\%`$ smaller than the VMC value. The error in the LDA energy density, which is obtained by weighting the energy of the LDA exchange-correlation hole by the density, is thus fairly insignificant at this point;  in contrast the LDA energy density in the atom core suffers from considerably larger deviations with the VMC.

Figure 1 also shows results for angle-averaged correlation and exchange holes at each point. In general, the exchange hole is the dominant contribution to the total hole. The correlation hole provides a correction to the exchange hole, reducing electron density near the reference electron and enhancing it towards the outside of the exchange hole, thereby making the total hole slightly deeper and more compact. At high densities this correlation hole is weaker and shorter-ranged than the LDA hole, a result in agreement with the prediction of gradient based models.  This trend however does not carry over to low densities, particularly in the atomic core \[Fig. 1(c)\]. The correlation hole here, while clearly shorter-ranged than the LDA hole, is also significantly larger and makes a much larger relative contribution to $`ϵ_{xc}(𝐫)`$. The WDA and ADA both do moderately well in predicting the magnitude and length scale of the correlation hole in this case.

In general the nonlocal density averaging methods fair less well at predicting correlation or exchange alone than the combination of the two. For the WDA in particular, the sum rules of the correlation hole and the exchange hole alone are not satisfied by satisfying the sum rule for exchange-correlation \[Eq. (14)\], with the correlation hole influenced more by the short-range behavior of the angle-averaged density and exchange the long-range behavior. This may contribute to the larger errors seen in $`n_c`$ and $`n_x`$ than for $`n_{xc}`$ observable in the cases with significant density variation.

Some sense of the general trends resulting from the different strategies for choosing the pair correlation function can be obtained from plotting the solution for the average interelectron distance $`\overline{r}_s`$ used by the nonlocal methods versus the $`r_s`$ obtained from the local density. The result for sample points at various densities throughout the unit cell is shown in Fig. 2. In the case of the ADA, for $`r_s`$ lower than the unit-cell average of 2.01 a.u., $`\overline{r}_s`$ hews closely to the local density value. It deviates from it at lower densities, either dramatically at the atom center (the one point below the average $`r_s`$ curve), or more gradually as one heads towards the interstitial point (the other low-density points plotted). The WDA values in contrast are all grouped about the unit-cell average value, regardless of location and density. In other words, the variation in the WDA exchange-correlation hole is primarily derived from variation in the nonlocal angle-averaged density, and very little from the variation of $`\overline{g}`$. These results reflect the character of the criteria used to determine the average density in each case. The weighting function of the ADA is an oscillatory function in real space, with a peak contribution at roughly 0.5 $`r_s`$ from the reference electron, so that a significant effect is observed only where the density variation about the electron is sufficiently rapid. In contrast, the weight implied by the WDA sum rule condition, $`u^2[1g^{heg}(u,r_s)]`$, is peaked near $`r_s`$ and has significant contributions out to 2 $`r_s`$, a distribution broad enough to to pick up the average density at almost every point in the crystal.

## V Exchange and correlation holes

As discussed earlier, the trends in the comparison of the DFT and VMC models for $`e_{xc}(𝐫)`$ are complex and difficult to characterize. Some insight into what is going on can be obtained by decomposing the exchange-correlation hole into exchange and correlation components and particularly by considering the response not only as a function of distance from the reference electron but including the complete information, $`n_{xc}(𝐫,𝐫+𝐮)`$, of the second electron’s position in the crystal.

The exchange hole is of use particularly in that it can be derived exactly and thus provide an unambiguous test of density functionals. Further, it is often the dominant part of the total exchange-correlation hole. It may be formally obtained as the exchange-correlation hole associated with the Slater determinant ground-state wave function of the noninteracting ($`\lambda =0`$) limit of Eq. 1. By applying this wave function in Eq. LABEL:eqnxc one has

$$n_x(𝐫,𝐫+𝐮)=\frac{1}{n(𝐫)}\underset{\sigma }{}\left|\underset{\alpha }{\overset{N_\sigma }{}}\psi _\alpha (𝐫)\psi _\alpha ^{}(𝐫+𝐮)\right|^2,$$
(20)
where $`\sigma `$ is a spin index and $`N_\sigma `$ the number of electrons with spin $`\sigma `$. In practice $`n_x`$ is obtained from LDA valence orbitals periodic on the $`3\times 3\times 3`$ simulation cell used in our VMC calculation. These produce a density that is indistinguishable from the VMC density within the statistical limitations of the Monte Carlo sampling. At $`𝐮=0`$ the exchange hole reduces to $`n(𝐫)/2`$, which reflects the Pauli exclusion prohibiting the occupation of the same point in space by two particles of the same spin. In addition, the exchange hole is strictly negative, as may easily be seen in the case of Si where the orbitals may be taken as real.

In Fig. 3 we show results for the exchange hole and angle-averaged hole weighted by distance for three positions near the bond center of the Si crystal. In the bottom half of the figure is a contour plot of $`n_x(𝐫,𝐫+𝐮)`$: the exchange hole at $`𝐫+𝐮`$ given a reference electron at the position $`𝐫`$. The plots cut through a chain of Silicon atoms in the (110) plane of the crystal; the reference electron’s position is shown as a white dot. The shading represents a change in value ranging from just less than zero for the white region to just above the absolute minimum of minus one-half the peak density of 0.087 a.u. in the black regions. The thick solid curve in the upper plots show the result of spherically averaging the exchange hole to obtain a function of distance from the reference point. The same weighting factor of $`2\pi u`$ as that of Fig. 1 is used. This result is compared to the various DFT models discussed previously.

In the case of the reference electron at bond center, the exchange hole is focused in the bonding region immediately surrounding the reference electron, with the bulk of the hole contained in the region between the atoms on either side of the bond, and a negligible weight on the nearest neighbor bonds. The electron density has been reduced by roughly 45$`\%`$ within the bond, i.e., in a diamond-shaped region extending 1.7 a.u. along the bond axis in either direction and 1.2 a.u. perpendicular to it. This accounts for 90$`\%`$ of the density of electrons with the same spin as the reference one. As one moves the reference electron either along the bond axis (c) or perpendicular (a) to it, the hole tends to remain focused on the nearby bond center, despite a small shift in weight towards the reference point. As before, more than 90$`\%`$ of the same-spin density is removed in the bond. This stiffness or insensitivity to the position of the reference electron is gradually lost as the reference electron is moved, say, to the atom center or towards the antibond point equidistant from three bonds. In the latter case, the exchange hole is again centered on the reference point and resembles a $`sp^3`$-hybrid atomic orbital.

The effect of the stiffness of the exchange hole on its angle average can be seen in the difference between the LDA and exact cases (thin and thick solid lines). In the bond center, Fig. 3(b), the LDA result closely matches the exact value. At position (a) off the bond axis, the bond center is oriented tangentially with respect to an angle average about the reference point and the exchange hole has become deeper and narrower than that obtained by the LDA. In contrast, in case (c), where the bond center is oriented radially out from the reference electron, and the electron is near to an atomic core, the hole has become shallower and wider than in the LDA. Although the two positions have similar densities and density gradients, the resulting deviations from the LDA are qualitatively different and point to a genuinely nonlocal behavior in the exchange hole. Nonetheless, neither the WDA nor the ADA are particularly sensitive to the position of the bond center with respect to the reference electron and provide no systematic improvement (or disimprovement) over the LDA fit.

The correlation holes corresponding to each case in Fig. 3 are shown in Fig. 4. Note the order of magnitude smaller range of density changes involved. The correlation hole $`n_c`$ takes on both negative values in regions where electrons are repulsed by the Coulomb potential of the reference electron and positive values representing regions where electron density has increased. The net electron density change is zero.

At the bond center, Fig. 4(b), the correlation hole contour plot can be roughly separated into two regimes: a deep and narrow well in the vicinity of the reference electron and a shallow longer-ranged response that may be described as a polarization of the nearest neighbor bonds. In the region of the pseudopotential core the correlation response is suppressed. The resulting angle-averaged hole has a complex and unusual structure as a function of distance. In comparison, the LDA and other DFT holes show a relative lack of structure and greatly overestimate the magnitude of the correlation energy per particle at bond center.

As the position of the reference electron is moved either off the bond axis (a) or along the axis (c), the correlation hole demonstrates a marked sensitivity to its position relative to the bond center. In either case the minimum is slightly offset in the direction of steepest increase in density. This is consistent with $`g_c`$ being nearly isotropic near the electron as expected from the cusp condition. Outside this short-range region the hole undergoes strong distortions in shape from that of the hole about the bond center. These longer wavelength features, as in the bond-center case, are well characterized in terms of the polarization of the nearby bonds. Bond polarization is particularly noticeable for the electron fixed on the bond axis (c) where there is a large shift of the electron density from the atom nearest the electron to the other side of the bond, with a peak on the bond axis opposite the reference electron. Nearest neighbor bonds are polarized in a similar fashion. In (a), the nearby bonds are polarized with density shifted from the nearer side of the bond with respect to the electron position to the farther, with the node of the polarization approximately along the bond axis. Consideration of the corresponding pair correlation function shows that these features in the correlation hole constitute a significant departure from an isotropic form of pair correlation function such as is assumed by the DFT models discussed in this paper.

Despite the reduction of information engendered by averaging over angle, the dramatic changes in the shape of $`n_c`$ with position show up in the angle-averaged hole, giving rise to nonlocal corrections to the LDA hole. The LDA model depicts a trend in $`n_c`$ to a broader and shallower hole as the density decreases from case (b) to (c) to (a) in the upper half of Fig. 4. The VMC hole in case (a) retains the character of the bond center hole, being narrower and weaker than that predicted by the LDA. The overall disagreement with LDA is less pronounced and some of the detailed structure is lost. As one moves along the bond axis to case (c), however, a significantly different trend occurs: the hole in the region of density reduction becomes deeper and broader relative to the LDA, and in particular now matches the width of the LDA hole. In addition to the stronger response in this region, a large peak appears at the hole edge in keeping with the zero sum rule of the correlation hole. This peak, at about 3.5 a.u., correlates with the position of electron density peaks on the nearest neighbor bonds in the contour plot. As with the case of exchange, these trends in $`n_c`$ poorly correlate to the changes in density or density gradient, and show a sensitivity to the larger nonlocal crystal environment.

Though the exchange and correlation holes both show significant and nonlocal deviations from the isotropic form typical of the homogeneous electron gas, these deviations are highly correlated with each other and thus tend to cancel in their sum. As the exchange hole tends to stay on the bond center, the bond polarization in the correlation hole is oriented so as to shift the center of the exchange-correlation hole towards the reference electron. A simple picture of this feature in the correlation hole is that it represents a response to the anisotropic distribution of the exchange hole by which some of the deviations of the exchange hole from an isotropic form are screened out in the correlation hole. This behavior may help to explain the relative success of the LDA in describing the exchange-correlation hole.

Neither the LDA nor the nonlocal DFT’s model with much fidelity the variation in the VMC correlation hole, aside from the general broadening of the hole with lower density. In particular the corrections of the WDA lead to a worse fit of both the correlation hole and its energy. The differences in exchange and correlation between WDA and VMC do tend to balance each other but not consistently: for example they cancel nicely in (a) and (b) but add in case (c). The net effect of these trends on the exchange-correlation energy is shown in Table II, where the total $`E_{xc}`$, $`E_x`$ and $`E_c`$ for several methods are compared. In addition to the methods discussed in detail in this paper, we present results from the PW91 version of the GGA  and those of a diffusion Monte Carlo (DMC) calculation.  The DMC method removes nearly all of the variational bias in the VMC correlation energy at a considerably larger computational cost. With the exception of the GGA, none of the density functional methods obtain a good estimate for $`E_x`$ or $`E_c`$ as compared to DMC, with errors of roughly 1 to 2 eV per atom. Estimates for $`E_{xc}`$ are much closer, especially for the WDA.

## VI Discussion

### A Limits of density averaging

The clear success of density averaging occurs at low densities where ADA and WDA obtain excellent exchange-correlation energy densities. WDA in particular predicts the shape of the exchange-correlation hole with exceptional fidelity, even in the extreme situation inside the atomic core. At high density, the small level of variation in the nonlocal angle-averaged density \[Eq. (18)\] limits the effects of density averaging to subtle alterations of the hole which do not provide a systematic improvement with respect to the LDA. More importantly, the discrepancy between the density-averaged holes and the VMC hole (much less the exact hole) are difficult to correlate with any known quantity.

However, density averaging must lead to a global and systematic improvement over the unit cell or else the quality of the total result may be disappointing. This is particularly the case of the ADA, which returns values very close to the VMC values at low density and slightly underestimates the magnitude of $`E_{xc}`$ at high density. Unfortunately, once the variational bias of the VMC energy is removed by a DMC calculation as shown in Table II this result proves to be a significant disimprovement with respect to the LDA. Moreover, although the different prescriptions for the WDA and LDA holes lead to extreme differences in the degree of nonlocality they incorporate, both end up with quite similar, reasonable predictions for the total $`E_{xc}`$. Apparently, the large error at low densities of replacing $`n(𝐫+𝐮)`$ with the local density $`n(𝐫)`$ in the definition of the LDA exchange-correlation hole is offset by correspondingly large changes in $`g`$, which are suppressed by density averaging in the WDA.

It is interesting to compare the density averaging approach to the gradient expansion used as the basis for the GGA. A notable result of the PW91 version of the GGA is that it not only provides a significant improvement of the LDA exchange-correlation energy, but of the exchange and correlation energies as well (Table II). Thus, it alone among the functionals we have studied shows promise to be a reasonable candidate as a correlation-only density functional for this material. In comparison, the WDA, although it returns a very good $`E_{xc}`$, has the worst estimate for $`E_c`$. However, the local density and its Laplacian, the formal equivalent of the information about the surroundings of an electron used in obtaining a GGA hole, constitutes but a small part of what is contained in the angle-averaged density $`n(𝐫,u)`$. The limitation of the WDA is the restriction to a pair correlation function form with only one variable parameter, one to which the hole in our case proves largely insensitive. It would be interesting if more of the information contained in $`n(𝐫,u)`$ could be used as input into a more flexible model for $`g`$, particularly for correlation or exchange separately. Possibly, the rate at which the density changes over the length scale of the hole, a factor which, as illustrated in Fig. 1, is clearly important in determining the overall shape of the WDA hole, could be used to influence the form of $`g`$ in analogy to the GGA.

Furthermore, $`n(𝐫,u)`$ has truly nonlocal information not accessible to GGA’s. A salient feature of the system in our study is the long-range limit of the angle-averaged density, which tends to a finite constant for a crystal and to zero for an atom or molecule. It is thus possible that low density points in these two systems have similar local or semilocal environments (and thus the same GGA holes) but significantly different exchange-correlation holes because of different long-range boundary conditions. Because of the nonlocal character of the WDA, these boundary conditions have a strong influence on the WDA hole, causing the length scale $`\overline{r}_s`$ to be nearly constant for the crystal but for the atom, to range from a finite value near the valence density maximum to infinity as the reference electron is removed from the atom. The WDA produces a reasonable fit to $`n_{xc}`$ at low density for both the Si atom and Si crystal, despite the markedly different $`\overline{r}_s`$ factors necessary to satisfy the particle sum rule.  It would seem, then, that the WDA should be capable to match or surpass the quality of the GGA – perhaps given a more accurate and flexible model pair correlation function form.

### B The exchange-correlation hole and energy density in the GGA

As discussed in Sec. II, the GGA models gradient expansion corrections to the LDA in terms of the gradient of the density alone, mapping corrections proportional to its Laplacian to a gradient squared term by an integration by parts in Eq. 7. This transformation has significant implications for the GGA models of $`n_{xc}`$. For example, each position in Fig. 1 is a critical point in the density, and as such characterizes the nature of $`n_{xc}`$ in its vicinity. At these points, the PBE model for $`n_{xc}`$ is indistinguishable from that of the LDA, a particularly serious error in the atom core region 1(c), where effects of inhomogeneity are most apparent. At the same time, corrections to $`n_{xc}`$ from the atom core and other regions where the Laplacian is large are mapped to the $`n_{xc}`$ of regions where the gradient is large, producing misleading corrections in these regions as well. The GGA is rather designed to provide gradient corrections to the LDA on a system-averaged basis. It does in fact capture general trends in exchange and correlation as discussed briefly in the next section, but clearly not on a local or point-by-point basis.

This approach is frequently justified by the observation that the energy density is not in itself a necessary component of an accurate DFT, rather only the total $`E_{xc}`$ and its functional derivative with respect to the density, the exchange-correlation potential.  Nonetheless, it is worth mentioning that the deviation of the LDA energy density from the VMC energy density closely follows the Laplacian of the density. For example, the minimum of the Laplacian, an hourglass-shaped region near the bond center, closely correlates with the region where we have found the largest positive deviation of the LDA energy density from the VMC value.  Likewise a large negative deviation in the LDA energy density occurs in the atom center where the Laplacian has a maximum, and a smaller negative deviation in the interstitial region where the Laplacian has a weak local maximum. This latter feature is clearly visible in the $`n_{xc}`$ and $`ϵ_{xc}(𝐫)`$ data of Sec. IV. We expect that a GGA model that fits both gradient and Laplacian terms could prove very useful in improving the energy density and thus also the exchange-correlation potential.

### C Orbital effects in exchange and correlation

Insight into the comparison of our VMC results and those obtained from the various models derived from the homogeneous electron gas may be obtained from a consideration of the semiconductor environment.

First of all, a distinguishing feature of systems such as molecules or semiconductors which have a finite energy gap is that the ground state may be described in terms of an exponentially decaying localized basis. In a periodic system these are the Wannier functions, which in Si can be defined by a unitary transformation from the four valence bands to four orbitals per unit cell localized about the bond centers $`𝐫_I`$

$$W_I(𝐫𝐫_I)=\underset{nk}{}U_{nk,I}\psi _{nk}(𝐫).$$
(21)
The symmetry of the crystal requires that each orbital $`W_I`$ must be related to the others by a space group operation. Such orbitals are well localized on a given bond and describe well the bonding character of the valence electrons. The exchange hole is most clearly represented in terms of these local functions and deviates strongly from the homogeneous electron gas model when these functions have little overlap. The exchange hole in terms of a Wannier basis is

$`n_x`$ $`(𝐫,𝐫+𝐮)=`$ (23)
$`{\displaystyle \frac{1}{n(𝐫)}}{\displaystyle \underset{\sigma }{}}\left|{\displaystyle \underset{I}{\overset{N_\sigma }{}}}W_I(𝐫𝐫_I)W_I^{}(𝐫+𝐮𝐫_I)\right|^2.`$
Using Wannier functions obtained from our LDA orbitals and the method of Marzari and Vanderbilt,  we find that 97$`\%`$ of the the on-top value $`n_x(𝐫_I,𝐫_I)`$ of the exchange hole at bond center is determined solely from the Wannier function localized on that bond site. Thus, the lack of sensitivity of the exchange hole with respect to reference position seen in Fig. 3 near the bond center is likely a reflection of the domination of the exchange hole in this region by a single Wannier function. This situation is reminiscent of that of the $`\mathrm{H}_2`$ molecule in which the exchange hole is constructed from a single orbital and is totally insensitive to electron position.

An exponentially decaying exchange hole typical of an insulator or finite system has the effect of lowering the exchange energy with respect to that of the LDA, owing to the more localized form of the hole as compared to that of the homogeneous electron gas; imposing a finite-ranged hole has been an ingredient in constructing successful GGA methods.  In our case the actual $`E_x`$ is 1.5 eV lower than the LDA result, in good agreement with the PW91 GGA result; likewise, the the exchange energy per particle $`ϵ_x(𝐫)`$ is lower than the LDA prediction through much, but not all, of the unit cell as shown in Fig. 5.

The semiconductor environment also affects the correlation hole in several ways. First, there is a finite energy cost to correlate electrons as compared to the homogeneous electron gas and therefore there should be a smaller correlation response and correlation energy than predicted in the LDA. This is seen in the VMC correlation energy density and energy per particle at almost all points in the unit cell, as shown in Fig. 5. Again this trend is consistent with the assumptions made in the PW91 version of the GGA but does not correlate with those made in the density averaging methods.

Secondly, as the electronic structure of the covalent bond plays an important role in exchange (in the form of a bond-orbital-like exchange hole near the bond center), it should do so as well in correlation. The signature of the existence of the covalent bond in the correlation hole is bond polarization, a clear feature of our VMC data. In the simplest bond-orbital picture,  this arises from the introduction into the noninteracting ground state of excited states in which a pair of electrons are excited from bonding to antibonding orbitals on the same bond site.  These states contribute a “left-right” correlation similar to that of $`\mathrm{H}_2`$ and other diatomic molecules in which the probability density of the second electron is shifted along the bond axis to the opposite side of the bond as that occupied by the first electron. A similar excitation of nearest neighbor electrons causes van der Waals like correlation between bonds.

These effects are evident in the “bond-right” case (c) of Fig. 4 where the antibond orbital has a peak. At the same time, the correlation hole associated with this left-right polarization is negligible when the reference electron is placed on the nodal plane of the antibond which passes through the bond center normal to the bond axis: that is, the situation of cases (a) and (b) of Fig. 4. Consequently, the weakness of $`n_c`$ relative to the homogenous electron gas based DFT models in these two cases and the enhancement of the correlation hole on the bond axis have a plausible explanation in the differing contribution of the antibond correlation to the holes in each case.  Moreover, bond polarization contributes noticeably to the functional variation of $`ϵ_c(𝐫)`$ in this bonding region. We find that the $`ϵ_c(𝐫)`$ obtained from VMC data has a maximum (that is, a minimum in the energy reduction obtained by electron correlations) on a narrow ridge centered on the antibond node passing through the center of the bond with normal parallel to the bond axis. Although some of the variation in $`ϵ_c(𝐫)`$ in the bonding region can be accounted for within the LDA, the behavior is sufficiently dissimilar to the local density variation as to lead to the largest discrepancy in the unit cell between the VMC and LDA results, as shown in Fig. 5.

The significance of the antibond correlation and bond polarization in general is that it provides a partial explanation of the cancellation of errors in the LDA exchange and correlation hole seen in their sum. The exchange-correlation hole of the homogeneous electron gas is essentially “dynamic” in nature, that is, responding to the position of the reference electron and largely insensitive to the details of electronic structure. The “nondynamic” or orbital-dependent features in the exchange hole, such as its insensitivity to particle position near bond center are a large potential source of error for the LDA and other models based on the homogeneous electron gas. The partial cancellation of these effects by a corresponding nondynamic feature of the correlation hole should then provide a combined hole much more amenable to the LDA and related density gradient corrections.

We have made preliminary estimates of the antibond contribution to the correlation hole $`n_c`$ and the correlation energy per particle $`ϵ_c(𝐫)`$ in the bonding region of Si using perturbation theory and perturbative configuration interaction methods  and localized orbitals derived from pseudopotential plane-wave DFT orbitals. We reproduce the qualitative differences in $`ϵ_c(𝐫)`$ between the points studied in this paper, as well as a reasonable reproduction of the bond polarization features of Fig. 4. Details of this calculation are to be reported in a later paper.

### D Error analysis

As discussed in Section III, the exchange-correlation hole suffers from fairly large plane-wave cutoff errors in the atomic core where the pseudopotential orbitals vary rapidly. In order to produce a reasonable correlation hole, particularly in the core, it was necessary to correct for this error by taking the $`n_x`$ and $`n_{xc}`$ expanded to the same cutoff. Then a best estimate for $`n_{xc}`$ was taken by adding the finite cutoff estimate of $`n_c`$ to the exact $`n_x`$, that is, an estimate of $`n_{xc}`$ obtained from correcting the exchange component of the Monte Carlo estimate:

$$n_{xc}=n_{xc}^{approx}+n_xn_x^{approx}.$$
(24)
The change in $`n_{xc}`$ is noticeable particularly for the atom center case Fig. 1(c), where the plane-wave correction results in a deeper and more localized $`n_{xc}`$. The resulting $`ϵ_{xc}(𝐫)`$, shown in Table I, varies from the raw VMC data by roughly $`8\%`$, increasing the agreement with WDA and ADA at the expense of the LDA value. However, given the poor agreement of the LDA and VMC values, the uncertainty in the VMC value does not significantly alter the assessment of a dramatic improvement by ADA and WDA in this region. A rather smaller change was noticeable for the angle-averaged hole for the “bond-right” position discussed in Figs 3 and 4. In this case, the difference is noticeable on the scale of the energy differences between the various theories, with the plane-wave correction favoring the ADA case. The correction was essentially negligible for the other points studied.

A second source of error in the VMC data is due to the statistical nature of the Monte Carlo estimates of expectations. Statistical sampling leads to roughly equal errors in the expectations of each measured plane-wave component of the pair correlation function, and a homogenous background noise (up to the resolution of the plane-wave expansion) in the real space behavior of the function. With the assumption of homogeneous and uncorrelated background noise, the weighted angle-averaged holes shown should have an error that is roughly constant at large interparticle distances. We find this to be a reasonable description of the long-range fluctuations in this quantity, resulting in a very small error for all the cases we have studied. On the other hand, the energy and particle sum rule of the exchange-correlation hole, both integrals over a large volume, suffer more serious cumulative effects from background noise, particularly at large interparticle distances where $`n_{xc}`$ is essentially zero. As a result, the VMC particle sum rule as estimated from the plots in Fig. 1 typically varies by 1 to 3$`\%`$ from the correct value. The values for the exchange-correlation energy per particle and energy density reported in our previous studies  have been made by adjusting the exchange-correlation hole at each $`\lambda `$ to obtain the correct particle sum rule, and numerical values at the various points considered in detail here are shown in Table I. This correction turns out to be noticeable at high density relative to the small differences between VMC and the various DFT models.

A final source of error is variational bias in the trial ground-state wavefunction resulting from its limited variational flexibility. The resulting discrepancy from the true ground state can be significant for correlation: the optimized VMC correlation energy is 14$`\%`$ higher than that estimated by the nearly exact DMC calculation. Exchange, on the other hand, is relatively unaffected by variational bias, as the VMC density and the associated single-particle orbitals are expected to be very accurate.

Convergence studies  of a Slater-Jastrow trial wavefunction for the Si atom reveal trends that we also expect to be observed in the present study. The correlation hole averaged over the position of the reference electron should be quite close to that of the true ground state, with deeper minima and sharper maximum in proportion to the change in $`E_c`$. The hole as a function of position in the system will in addition show subtle alterations in shape, while preserving basic qualitative features, for example the bond polarization in Fig. 4.

The correlation hole can be expected to be more accurate at higher densities which carry greater weight in a variational optimization of the wavefunction; however the Slater-Jastrow wavefunction guarantees important conditions on the hole at any density, such as the cusp condition and sum rules, so a reasonable estimate of the hole should be obtained even at quite low densities. Assuming that the generic trend at high density is to deepen the correlation energy of the hole by roughly $`14\%`$, we expect the overestimate of $`ϵ_c(𝐫)`$ by the DFT models, as shown for example in Fig. 5 to be reduced by between 10 and 50$`\%`$ in the Si bond. Qualitative trends in this region will not be changed if the error in $`ϵ_c(𝐫)`$ does not vary dramatically with position.

In general, the various errors in the VMC data (plane wave cutoffs, sum rule enforcement, and variational bias) are most significant for the angle-averaged $`n_{xc}`$ and exchange-correlaction energy per particle at high density. Even though estimated corrections are rather small, the agreement with the various DFT’s studied is quite close and even small changes in the VMC data can be significant. Consequently, the errors in our data contribute to the general difficulty we find in assessing the high density trends of the WDA and ADA. On the other hand, the difference between the DFT theories and the VMC are larger for the angle-averaged correlation hole. The qualitative trends in $`n_c`$ near bond center discussed in this paper are accordingly unaffected by error corrections though there are small differences in the quantitative value of $`ϵ_c(𝐫)`$ at the “bond-right” position. Finally, the improvement of ADA and WDA over LDA at low density is so dramatic that the observed error corrections to the VMC are small in comparison.

## VII Conclusion

We have carried out a detailed analysis of the exchange-correlation hole and energy per particle in the Si crystal, comparing various DFT models to accurate numerical data calculated with the VMC method.

We find that the WDA and ADA help overcome the major defects of the local density approximation at low densities and especially in the pseudopotential atom core, where the rapid change in density relative to the length scale of the LDA hole dramatically affects the shape and range of the hole, and thus improve the fit with the VMC energy density. The remaining discrepancies are due to the failure of density averaging to provide significant information at intermediate and high densities, where the inhomogeneity in the density is less severe, but the contribution to the energy density of subtle effects of the inhomogeneity is noticeable because of the higher density. These discrepancies are, at least in principle, the result of the inflexibility of the scaling form of the pair correlation function used to fit the hole, which is insensitive to subtle variations in density.

The detailed investigation of exchange and correlation at high density reveal the importance of orbital correlations in both cases. We find that the exchange hole is well described by a Wannier bonding state near the bond center, and that a “left-right” correlation or bond polarization plays an important role in the spherically-averaged correlation hole and has a noticeable effect on the correlation energy density. This bond polarization has the effect of countering the coarse-scaled deviation of the exchange hole from an isotropic form centered on the electron. As a result, the exchange-correlation hole and energy density is much more reliably fit by the DFT models studied than either exchange or correlation alone.

Viewed in another sense, these orbital effects lead to serious defects in the correlation hole and energy of models derived from the homogeneous electron gas. The WDA and ADA conspicuously fail to improve on the LDA both in the system averaged energies $`E_x`$ and $`E_c`$ and in detail, particularly in the bonding region. Here, the nonlocal corrections to the LDA employed by these models fail to account for the physical features responsible for the most significant errors in the LDA model.

Of all the methods considered in this paper, the GGA alone accurately predicts $`E_x`$ and $`E_c`$ for Si. Despite the lack of a direct comparison between our data and current GGA models for $`n_{xc}`$ the assumptions made in constructing the GGA exchange and correlation hole appear to be generally born out by our results, with a notable exception in the low density, large inhomogeneity limit for the correlation hole. Thus we expect that the accurate GGA results for $`E_x`$ and $`E_c`$ are due to an improved description of $`n_x`$ and $`n_c`$ on average, if not locally. Interestingly, the point by point trends of our data seem to be best characterized by the Laplacian of the density, in contrast to the approach taken by current gradient-only GGA models. Moreover, exchange and correlation taken alone show nonlocal qualitative features at high density that may prove difficult to describe in terms of any local density expansion.

It would be of great value with respect to augmenting exact-exchange calculations to have a compact expression of the bond polarization features of the correlation hole in terms of localized atomic or bond-centered orbitals, possibly with a short-range correction in LDA or GGA. We have done preliminary calculations on modeling this bond polarization in terms of excitations of electrons into antibonding states. However a conveniently usable and compact form for the correlation hole at long range remains to be found.

This work was supported in part by the Department of Energy (Grant No. DE-FG02-97ER45632) and the National Science Foundation (Grant No. DMR-9724303). It was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48. We wish to thank Nicola Marzari for his kind help in providing us Si Wannier orbitals.