Energy density in density functional theory: Application to crystalline defects and surfaces

# Energy density in density functional theory: Application to crystalline defects and surfaces

Min Yu Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801    Dallas R. Trinkle Department of Materials Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801    Richard M. Martin Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801
July 19, 2019
###### Abstract

We propose a method to decompose the total energy of a supercell containing defects into contributions of individual atoms, using the energy density formalism within density functional theory. The spatial energy density is unique up to a gauge transformation, and we show that unique atomic energies can be calculated by integrating over Bader and charge-neutral volumes for each atom. Numerically, we implement the energy density method in the framework of the Vienna ab initio simulation package (vasp) for both norm-conserving and ultrasoft pseudopotentials and the projector augmented wave method, and use a weighted integration algorithm to integrate the volumes. The surface energies and point defect energies can be calculated by integrating the energy density over the surface region and the defect region, respectively. We compute energies for several surfaces and defects: the surface energy of GaAs, the mono-vacancy formation energies of Si, the surface energy of Au, and the interstitial formation energy of O in the hexagonal close-packed Ti crystal. The surface and defect energies calculated using our method agree with size-converged calculations of the difference between the total energies of the system with and without the defect. Moreover, the convergence of the defect energies with size can be found from a single calculation.

## I Introduction

Total energy is one of the most important quantities in a solid state system, as it determines the stable configuration, and its derivatives provide other equilibrium properties. The assignment of energy to particular finite volumes can provide additional detailed information, such as defect formation energies. For example, the surface energy of a given facet of a crystal is meaningful for predicting the equilibrium crystal shape and preferred crystal growth directions, and should depend only upon the properties of the surface. The formation energy of a point defect is important to understanding the phase stability, and should depend only upon the properties in the vicinity of the defect. However, density functional theory calculationsKohn and Sham (1965) provide only a single total energy for a given configuration, not a spatial partitioning of energy into additive contributions—instead, defect energies are determined as a difference of two of more separate total energy calculations.

The energy density methodChetty and Martin (1992) can provide the formation energies for more than one point defect, surface or interface in a single calculation, as well as a picture of the distribution of the energy among the surrounding atoms. The energy density formula derived by Chetty and Martin is a reciprocal-space expression for norm-conserving pseudopotentialTroullier and Martins (1991) (NCPPs) with local density approximationCeperley and Alder (1980); Perdew and Zunger (1981) (LDA), where the ion-ion interaction is treated in a similar manner as the Ewald sum.Ewald (1921) The energy density is not unique since there are multiple definitions of the kinetic and Coulomb energy densities that integrate to the same well-defined total energy. The different definitions can be considered as gauge variability; defining a unique gauge-independent energy requires the identification of spatial volumes where the gauge differences integrate to zero. Chetty and MartinChetty and Martin (1992, 1992) showed that the surface energy for a crystal can be calculated by an integral over a region to high-symmetry planes within the bulk of the crystal, where symmetry ensures that the gauge dependent terms integrate to zero. Therefore, two polar surface energies such as the and the surfaces of a zincblende semiconductor GaAs can be integrated independently in a single calculation. This compares with more complex approaches that use multiple wedge-geometries to extract polar surface energies for specific geometries.Zhang and Wei (2004) Rapcewicz et al.Rapcewicz et al. (1998) followed the energy density formalism, but generalized the method to low-symmetry system such as the surface of GaN and the surface of SiC by introducing Voronoi polyhedra for the integration volumes; however, it should be noted that Voronoi polyhedra are gauge-independent integration volumes only in specific situations. RamprasadRamprasad (2002) extended the application of energy density method from surfaces to point defects of metals and presented two applications on the monovacancy of Al and the surface of Al.

We reformulate the expression of the real-space energy density for the projector-augmented waveBlöchl (1994) (PAW) method and the norm-conserving and ultrasoft pseudopotentialsVanderbilt (1990) (USPP), grouping terms in a similar manner. We implement the energy density method in the Vienna ab initio simulation packageKresse and Furthmüller (1996); Kresse and Joubert (1999) (vasp). We also demonstrate the usefulness of assigning an energy to each atom using gauge-independent integration volumes. Excluding cases that can be determined by symmetry—such as the energy of an atom in the diamond structure that is half the total energy of a unit cell—the assignment is not unique. Nevertheless, there are two primary reasons for developing an approach based on gauge-independent volumes. One advantage is that it provides an automatic procedure to choose volumes that do not depend upon the choices of planes or polyhedra in the calculations mentioned above. The resulting surface and defect energies is the same as for the special cases above, but can also be used for general cases. In addition, the energy per atom partitions space to individual defect regions, and provides a measure of finite-size errors in a single calculation. Moreover, within an individual defect—e.g., a vacancy in silicon—the atomic energy shows the contribution of individual atomic relaxations to the total formation energy, which can be useful for understanding changes to the stability of individual defects and surfaces.

The assignment of energy to an atom derives from Bader’s “atoms in molecules” theory.Bader (1990) Although our work is different, one element is the same: the calculation of volume around each atom where the kinetic energy is unique. We extend this to find a charge neutral volume for a unique classical Coulomb energy. In Bader’s work, a local form of virial theorem relates the local electron kinetic energy to the potential energy density, and defined atomic energies within each volume bounded by zero flux surface of the gradient of electron density. However, this results in charged units (so-called “Bader charges”) which have long range forces. Thus it is arbitrary to assign the energy to one region when it is in fact shared by interacting regions. Also, the Bader approach does not consider the exchange-correlation energy of density functional theory. A recent paperRodríguez et al. (2009) has reported a way to include such terms in a form of the virial theorem; they found different values from those given by the original method, making the applicability of a local virial theorem to density functional theory questionable. To overcome these difficulties, we instead define two integration volumes: one based on kinetic energy, and another based on potential (electrostatic) energy.

We derive the energy density methodology for PAW, and apply it to several defect types in solid state systems. Section II reviews the energy density expression derived by Chetty and Martin for NCPPs, and derives the reformulated expression for the PAW method. The energy density contains a gauge-dependent kinetic energy density, a gauge-dependent long-ranged classical Coulomb energy density, a well-defined exchange-correlation energy density, and short-ranged terms grouped in an on-site energy for each ion. We also consider the method to integrate the gauge-dependent terms. Section III presents four applications: the surface energy of GaAs, the mono-vacancy formation energy of Si, the convergence test of surface energy of Au, and the interstitial formation energy of O in the hexagonal close-packed Ti crystal. These examples show the generality of the method, and the new information extracted such as convergence of formation energy with size and the spatial partitioning of defect energy.

## Ii Methodology

In density functional theory, the standard form for the total energy of a crystal is usually written in reciprocal space. There are many forms convenient for the total energy written in terms of the wavefunctions and/or eigenvalues.Martin (2004) For our purposes, we consider the expressions in real space where the Kohn-Sham wavefunctions are used only for the kinetic energy and for non-local terms in the pseudopotential. In the norm-conserving pseudopotential approximation, the total energy is, in atomic units (),

 Etot=−12∑nkfnk∫dr~ψ∗nk(r)∇2~ψnk(r)+EXC[ρe(r)]+∑μEnlμ+∫drρe(r)∑μVlocμ(r−Rμ)+12∫drρe(r)VH(r)+∑μ<νZμZνRμν (1)

where and are the valence pseudo-wavefunction and the electron occupation number for the band, for wavevectors within the first Brillouin zone, and with valence electron density . The first term is the independent electron kinetic energy. The second is the exchange-correlation energy , where is the exchange-correlation energy per electron; in the local density (LDA) or a generalized gradient (GGA) approximation it is a function of the density or the density and its gradient. The fermion nature of many-body interacting electrons is approximated by this exchange-correlation potential. The third term is the energy due to the non-local part of the pseudopotential

 Enlμ=∑nk∑ℓ∫dr~ψ∗nk(r)Vnlμℓ(∣∣r−Rμ|)℘ℓ~ψnk(r), (2)

where is the th component of the non-local pseudopotential, with the projection operator on angular momentum . This term is site-localized (non-zero only within the core radius around a site) so that the total energy involves a sum over the sites at position . The last three terms of Eqn. (1) are the long-ranged Coulomb interactions. The fourth and fifth terms are the interaction of the electrons with the local ionic pseudopotential and with themselves that can be written as one-half the interaction with the Hartree potential . The sixth term is the valence charge of ion at position —the ion-ion Coulomb interaction energy that is the same as for point charges since the cores are assumed not to overlap. Finally, there can be non-linear core corrections not shown here, but which can be expressed in terms of involving the core density similar to the PAW method.

We use the PAW methodBlöchl (1994) for calculations, where the total energy has a similar form

 Etot=−12∑nkfnk∫dr~ψ∗nk(r)∇2~ψnk(r)+EXC[~ρ+^ρ+~ρc]+∑μ(E1μ−~E1μ)+EH[~ρ+^ρ]+∫VH[~ρZc](~ρ+^ρ)dr+∑μ<νZμZνRμν. (3)
 E1μ=∑ijρij⟨ϕi|−12∇2|ϕj⟩+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯EXC[ρ1+ρc]+¯¯¯¯¯¯¯¯¯¯¯¯¯EH[ρ1]+∫VH[ρZc](ρ1)dr~E1μ=∑ijρij⟨~ϕi|−12∇2|~ϕj⟩+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯EXC[~ρ1+^ρ+~ρc]+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯EH[~ρ1+^ρ]+∫VH[~ρZc](~ρ1+^ρ)dr (4)

for . Quantities with a tilde are obtained by pseudization, and a superscript 1 for quantities evaluated inside atom-centered spheres on a radial grid. For each atom-centered sphere, the pseudo-partial waves match all-electron partial waves at the sphere boundary and outside the augmentation region. The smooth projector functions are dual to the pseudo-partial waves, and are the occupancies of augmentation orbitals . Then is the soft pseudovalence electron density, and are the on-site charges (full and pseudized) localized around each atom, is the compensation charge, and are the frozen core charges (full and pseudized), and are the sum of the nuclei and core charges (full and pseudized). The electrostatic interactions—electron-electron, electron-ion, and ion-ion interactions (last three terms in Eqn. (3))—are collectively the “classical Coulomb” energy. The short-ranged terms for individual ions are . We derive the total energy density for the PAW and pseudopotential methods as

 e(r)=t(r)+eCC(r)+eXC(r)+Eon-siteδ(r−Rμ) (5)

and use gauge-independent integration over Bader and charge-neutral volumes to define atom-centered energies.

### ii.1 Kinetic energy density

The kinetic energy density is gauge dependent and can be expressed as asymmetric or symmetric functional,Martin (2004)

 t(a)(r)=−12∑nkfnk~ψ∗nk(r)∇2~ψnk(r)t(s)(r)=12∑nkfnk|∇~ψnk(r)|2. (6)

The difference between asymmetric and symmetric kinetic energy density is a gauge-dependent term proportional to the Laplacian of pseudo electron density,

 (7)

The integral of the two forms of kinetic energy density is equal when the gauge-dependent integral vanishes; e.g., for infinite or periodic systems. In Section II.5, we will integrate over a discrete set of atom-centered volumes where the gauge-dependent integrals also vanishes—hence, uniquely defined kinetic energies for atoms in a condensed system. Note that continuous wavefunctions can have cusps in their gradient, thus the asymmetric form of the kinetic energy density can be ill-defined. Chetty and Martin chose the symmetric form for the kinetic energy density as it appears in the variational derivation of the Schrödinger’s equation,Slater (1937) and hence is a more fundamental quantity. However, the kinetic energy density is unique except for terms proportional to the Laplacian of pseudo electron density; if we integrate over volumes where the gauge-dependent term of Eqn. (7) (c.f. Section II.5), then either form of Eqn. (6) gives the same kinetic energy. For a planewave basis, the asymmetric kinetic energy density is well-defined everywhere—i.e., there are no cusps in the wavefunction gradient—and is computationally less demanding to calculate. In the PAW method, the total kinetic energy density contains three terms

 t(a)(r)=~t(a)(r)+t1(a)(r)−~t1(a)(r). (8)

The first term, is also the first term in Eqn. (3) and is expressed by using the pseudo wavefunction. The on-site kinetic energies and are first terms in Eqn. (4), and are included in the short-ranged on-site energy .

### ii.2 Classical Coulomb energy density

The total classical Coulomb energy of a system with electrons and nuclei can be written as

 ECC=12∫dr∫dr′ρe(r)ρe(r′)|r−r′|+∫drρe(r)∑μVlocμ(r)+∑μ<νZμZνRμν (9)

where is the sum of soft pseudoelectron density and compensation charge , and are representing different nuclei and is the distance between two nuclei. There are various ways to calculate the electrostatic energy.Martin (2004) In the Ewald method,Ewald (1921) terms are grouped with a Gaussian charge density around each atom so that the sum can be calculated by sums in both real and reciprocal space. However, this is not useful for constructing an energy density in real space. Instead, methods that involve only smooth densities for each ion can be used to construct expressions for the Coulomb energy that are expressed only in real space.

#### ii.2.1 Smeared ions

We introduce a fictitious localized charge distribution , which gives rise to a local pseudopotential (c.f. Section F.3 of Martin, 2004) for ion as

 ρlocμ(r−Rμ)=−14π∇2Vlocμ(r−Rμ). (10)

The Coulomb interaction energy between two ions and is

 Elocμν(|Rμν|)=ZμZνRμν=∫drρlocμ(r−Rμ)Vlocν(r−Rν), (11)

and the self energy on each ion is

 Eselfμ=12∫drρlocμ(r)Vlocμ(r). (12)

The total classical Coulomb energy of a system with electrons and nuclei can be written as

 ECC=12∫dr∫dr′ρe(r)ρe(r′)|r−r′|+∫drρe(r)∑μVlocμ(r)+∑μ<νZμZνRμν=∫dr18π|∇Vtot% (r)|2−∑μEselfμ (13)

with total classical Coulomb potential . The Hartree, local, and ion-ion interaction terms (last three terms of Eqn. (3)) are combined into the classical Coulomb term.

We use the Maxwell energy for the total electrostatic component of the energy density. As a fictitious charge density defined in Eqn. (10), the total neutral charge density . With this definition, the Maxwell form of the classical Coulomb energy density can be written as

 eMaxwellCC(r)=18π|∇Vtot(r)|2. (14)

Similar to the kinetic energy density, the classical Coulomb energy density is unique up to a gauge transformation. The asymmetric form of the classical Coulomb energy density is

 e(a)CC(r)=−18πVtot(r)∇2Vtot(r)=12Vtot(r)ρtot(r), (15)

and the gauge-dependent term is the difference of Eqn. (14) and Eqn. (15),

 e(a)CC(r)−eMaxwell% CC(r)=−18π∇⋅[Vtot% (r)∇Vtot(r)]. (16)

As with the kinetic energy density, we can obtain a gauge-independent classical Coulomb energy as an integral over any volume bounded by a zero-flux surface of the gradient of the total Coulomb potential.

The PAW method introduces the soft compensation-charge , and the Hartree energy is

 EH=~EH+(E1H−~E1H)=EH[~ρ+^ρ]+∑μ¯¯¯¯EH[ρ1]−∑μ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯EH[~ρ1+^ρ] (17)

The first term is related to the soft valence-charge density and the soft compensation-charge density, and is included in the classical Coulomb energy density. The last two short-ranged terms are related to the short-ranged on-site energy , as are the electron-ion interactions.

#### ii.2.2 Model charge density

In practice, the local charge density due to local pseudopotential can vary rapidly (c.f. Figure 1), which causes numerical errors in a real calculation. To improve our numerical accuracy, we then introduce a model charge density to solve this problem. The model charge density is chosen as obeying following constraints: a spherically symmetric functional, which is centered at each ion, zero beyond the cutoff radius of local pseudopotential, and normalized as negative to the local charge density within the cutoff radius. The total charge density can be rewritten as

 ρtot(r)=ρloc(r)+ρmodel(r)+ρe(r)−ρmodel(r)=∑μ[ρlocμ(r)+ρmodelμ(r)]+δρ(r), (18)

where is a neutral and spherical charge density for each ion; is the difference between the valence electronic charge density and the model charge density. The asymmetric form of the classical Coulomb energy is

 ECC[ρtot]=ECC[ρloc+ρ% model]+∫dr(Vloc+Vmodel)δρ+12∫drδV% H[δρ]δρ−∑μEselfμ. (19)

The first term is

 ECC[ρloc+ρmodel]=∑μ<νEloc+modelμν(|Rμν|)+∑μEloc+modelμ. (20)

The electronic interaction between two neutral atoms is zero when there is no charge overlap, since all moments are zero for spherical charge distributions. Therefore, the first term in Eqn. (20) is zero. Combining the second term with the self energy in Eqn. (12), we have

 Eloc+modelμ−Eselfμ=12∫drVmodelμ(r)ρmodelμ(r)+∫drVlocμ(r)ρ% modelμ(r), (21)

which is a constant for each species of ions and can be canceled when studying the defect energies. Neglecting this constant term, the asymmetric form of the classical Coulomb energy density in Eqn. (19) is

 eCC(r)=[Vloc(r)+12VH(r)+12Vmodel(r)][ρe(r)−ρmodel(r)], (22)

where , and are already known in real space. Different model potential and model charge density can be constructed as long as they obey above constraints. In this work, the model charge density is a polynomial functional with continuous zeroth-, first- and second-order derivatives at and ; for ,

 ρmodelμ(r)={215πr3c[1−10u3+15u4−6u5]:rrc (23)

As shown in Figure 1, the local charge density varies rapidly with respect to radius , while the model charge density parametrized in Eqn. (23) smoothly decays to zero as increasing to . The corresponding potential is

 Vmodelμ(r)={15rc[12−14u2+28u5−30u6+9u7]:rrc (24)

The model charge density gives faster numerical convergence on a regular spatial grid compared to the rapidly varying local charge density. This model charge density has been tested by calculating the surface energy based on a Si bulk (8 atom) and a Si (111) slab (16 atom). The energy difference between the values calculated from total energy calculation and from energy density integration is about 7meV. It has also been tested on O atoms, and O molecules for various grid sizes, where we had a convergence problem by using a rapidly-varying local charge density. For a wide range of grid sizes, the total energy calculations converge to a precision of 1meV, while the difference between results calculated from those two methods is up to 0.4eV. However, the energy difference can be reduced below 1meV with the smooth model charge.

### ii.3 Exchange-correlation energy density

The exchange-correlation energy of many-body interacting electrons can be expressed in terms of an exchange-correlation hole which tends to be localized around each electron. In density functional theory it is usually treated as a function of the local density and its gradients, which is determined by the choice of exchange-correlation functional. For both the local density approximation (LDA) and generalized gradient approximation (GGA), the gauge-independent exchange-correlation energy density is

 eXC(r)=ρe(r)εXC[ρe(r),|∇ρe(r)|], (25)

where is exchange-correlation energy per electron.

### ii.4 On-site energy

The last term of the energy density in Eqn. (5) is short-ranged. For the PAW method, the on-site energy for each ion is composed of kinetic energy, exchange-correlation energy, and Coulomb energy including electron-electron and electron-ion interactions in the augmentation region, and is

 Eon-site=(E1μ−~E1μ)δ(r−Rμ), (26)

with the on-site energies and expressed in Eqn. (4). In practice, we calculate and for each ion using radial grids. For NCPPs and USPPs, this short-ranged term corresponds to the non-local pseudopotential energy. For USPPs,

 Enlμ=∑nk∫dr~ψ∗nk(r)(∑ijDionij|βi⟩⟨βj|)~ψnk(r), (27)

where the coefficients and the projector functions vary depending on the atomic species. For NCPPs,

 Enlμ=∑nk∑ℓ∫dr~ψ∗nk(r)Vnlμℓ(|r−Rμ|)℘ℓ~ψnk(r). (28)

### ii.5 Gauge dependence and uniqueness

We have two energy density terms to be integrated which are gauge dependent: the kinetic energy density and the classic Coulomb energy density. Defining a gauge independent energy requires integrating these energy densities over volumes to cancel out any gauge dependence. Previously, Chetty and MartinChetty and Martin (1992) integrated over Wigner-Seitz cells in a supercell; Rapcewicz et al.Rapcewicz et al. (1998) constructed a Voronoi polyhedron for each comprised atom. However, these volumes are not the general solution to removing gauge dependence. For the kinetic energy density, the gauge dependence, Eqn. (7), is proportional to the Laplacian of electronic charge density. Hence, we integrate over a volume where the gradients of electron density has zero component along surface normal direction , —the zero-flux “Bader” volume.Bader (1990) For the classical Coulomb energy density, the gauge dependence, Eqn. (16), is proportional to the Laplacian of the potential. Hence, we integrate over a volume where the electrostatic field has zero component along surface normal direction , —the zero-flux charge-neutral volume.

In this work, we construct two different volumes: the Bader volume is used to integrate kinetic energy density and exchange-correlation energy density, and the charge-neutral volume is used to integrate classical Coulomb energy density. Each of these volumes is “atom-centered”—that is, it contains one atom somewhere in the volume—and they each partition space: the union of all volumes is the total supercell volume, and the intersection of any two volumes is zero. We define these volumes on the same regular spatial grid used to represent the charge density and energy density terms. Accurate definition of the volumes and integration uses a weighted integration schemeYu and Trinkle (2010) that has quadratic convergence in the grid density. Finally, we use the integral over the gauge-dependent kinetic (Eqn. (7)) and classical Coulomb (Eqn. (16)) terms as an estimate of the integration error for atomic energies. This error estimate has a sign, so it is possible for the magnitude of the error of the energies of neighboring volumes and to be greater than the magnitude of error integrating over .

### ii.6 Summary

We summarize the procedure of calculating atomic energy using energy density method in Table 1. The total energy density of Eqn. (5) contains a gauge-dependent kinetic energy density, a gauge-dependent classical Coulomb energy density, a gauge-independent exchange-correlation energy density, and a short-ranged non-local energy. The well-defined atomic energy can be calculated with two different integration volumes. The Bader volume is employed for the integral of the kinetic energy and the exchange-correlation energy, and the charge neutral volume for the integral of the classical Coulomb energy.

## Iii Applications

To verify our implementation of the energy density method and highlight the new information it reveals, we perform DFT calculations with vasp on the GaAs(110) surface, Si monovacancy, Au(100) surface and O interstitial in Ti. We integrate the energy densities around the defect regions, and compare the integrated defect energies with values given by total energy calculations and experiments. Finally, the convergence of the atomic energy to bulk values within a single calculation shows the convergence (or lack of) for each calculation.

### iii.1 GaAs(110) surface

The GaAs(110) surface contains equal numbers of Ga and As atoms: a stoichiometric or non-polar surface. The surface energy of a stoichiometric slab is

 γsurf=12A(Eslab−NslabEbulkNbulk), (29)

for surface area , where is the total energy of a GaAs slab with pairs of GaAs atoms, and is the total energy of GaAs bulk with pairs of atoms. Our DFT calculations are performed with the PAW method,Blöchl (1994) with the local density approximation (LDA)Ceperley and Alder (1980); Perdew and Zunger (1981) for the exchange-correlation energy. The valence configurations for Ga is with cutoff radius 1.01Å, and As is with cutoff radius 1.11Å; this requires a plane-wave basis set with cutoff energy of 650eV. This gives a lattice constant of 5.6138Å for zincblende GaAs, compared with the experimental lattice constant of 5.65Å. The supercell contains 11 layers of atoms with a pair of GaAs atoms on each layer, and a vacuum gap of 8Å to prevent the interaction between slabs under periodic boundary conditions. We use Monkhorst-Pack k-point meshesMonkhorst and Pack (1976) of for bulk eight-atom cells, and for the slab supercell; Brillouin-zone integration uses Gaussian smearing with for electronic occupancies, and the total energy extrapolated to . We represent the charge density and compute energy densities on a grid of . Geometry is optimized to reduce forces below 5meV/Å. This gives a surface energy of 50meV/Å; this agrees with Moll et al.’s valueMoll et al. (1996) of 52meV/Å, Qian et al.’s valueQian et al. (1988) of 57meV/Å, Choudhury et al.’sChoudhury et al. (2008) LDA value of 50meV/Å, and the experimental valueMessmer and Bilello (1981) of meV/Å.

Figure 2 shows the energy change from bulk for each layer by integrating over volumes that eliminate gauge dependence. The change in energy shows differences from bulk that are mainly confined to the first two layers; the bulk-like response of the interior layers—not just for the total energy, but also the individual contributions to the energy. Determining the size-convergence of a surface calculation with total energy alone requires computing surface energies for multiple sizes; in our case, the bulk-like behavior of our center layers indicates a small finite-size error without requiring multiple size calculations. We can integrate the surface energy by adding the energies from the first two layers; our surface energy is meV/Å, which agrees well with the total-energy calculation of surface energy. The error estimate is specifically for the integration error over the Bader and charge-neutral volumes. Note also that we can compute the energy of each surface independently; for surfaces with different chemistry, this allows for two surface energies to be calculated from a single supercell.

Figure 3 shows the Bader and charge-neutral volumes for Ga and As atoms in a plane of GaAs. The Ga and As atoms all lay in a plane, and the intersection of the surfaces show the difference between the two atom-centered volumes. The Bader volumes have zero flux of the gradients of charge density through their surfaces, and are used to integrate a unique kinetic and exchange-correlation energy. The charge-neutral volumes have zero flux with the total electrostatic field, and are used to integrate a unique classical Coulomb energy. These volumes are different also from the Voronoi volumes around each atom. The atomic volumes, like the individual components of energy, become bulk-like in the center of slab. Atoms at the free surfaces have volumes that extend into the vacuum. Besides the different surfaces, the Bader volumes of As are larger than the As charge-neutral volumes.

### iii.2 Si monovacancy

The monovacancy in bulk Si is a simple point defect in a semiconductor, which has been studied theoretically and experimentally. From total energies, the formation energy of a vacancy is

 ΔHv=EN−1v−N−1NEN, (30)

where and are the total energy of the and atom supercells with and without one vacancy. WrightWright (2006) performed LDAPerdew and Zunger (1981) and GGA-PBEPerdew et al. (1996) calculations in 215, 511-, and 999-atom supercells to get formation energies of 3.53eV, 3.49eV, and 3.47eV with LDA and 3.66eV, 3.63eV, and 3.62eV with GGA. Puska et al.Puska et al. (1998) performed LDA calculations in 31-, 63-, 127-, and 215-atom supercells to get formation energies of 3.98eV, 3.42eV, 3.44eV, and 3.31eV. Experiments have found a formation energy of eV.Dannefaer et al. (1986)

Our DFT calculations are performed with the PAW method with the generalized gradient approximation (GGA) of Perdew and Wang (PW91)Perdew and Wang (1992) for the exchange-correlation energy. The valence configurations for Si is with cutoff radius 1.01Å; this requires a plane-wave basis set with cutoff energy of 417eV. This gives a lattice constant of 5.4674Å for diamond Si, compared with the experimental lattice constant of 5.43Å. The simple cubic supercell with a vacancy contains 63 atoms. We use a Monkhorst-Pack k-point mesh; Brillouin-zone integration uses Gaussian smearing with , and the total energy extrapolated to . We represent the charge density and compute energy densities on a grid of . Geometry is optimized to reduce forces below 5meV/Å. This gives a formation energy of 3.65eV.

Figure 4 shows the energy change from bulk for shells surrounding a Si vacancy. The primary contribution to the vacancy formation energy comes from the first five shells, becoming bulk-like at larger distances. Summing the atomic energies up to the fifth shell gives a formation energy of , which is similar to the total energy calculation. The importance of the fifth shell over the third and fourth shells can also be seen in charge disturbances from a vacancy. KaneKane (1985) showed charge disturbances around a Si monovacancy out to the 27 shell, while the first two shells contribute 60% of the charge disturbance. The most striking future was charge concentrates on the planar zigzag chains of atoms, such as [000], [111], [220], [331], [440], as so on. Twelve such chains exist by symmetry. After reaching the fifth shell at , the charge decays monotonically along the coplanar chains. Kane connected this result to the importance of the fifth-neighbor interaction in the valence force modelMcMurry et al. (1967) of covalent phonon spectra. The valence force model is an empirical model connecting force constants to the electronic configuration. The fifth-neighbor interaction is proportional to from changes in bond angles and along a zigzag chain. The fifth neighbor has a stronger interaction with the bond-bending than third and fourth neighbors; we see a similar change in energy for the vacancy.

### iii.3 Au(100) surface

Au(100) is a low-index metallic surface without a reconstruction. To calculate the surface energy with Eqn. (29), the thickness of the slab must be increased until the surface energy converges. Our DFT calculations are performed with the PAW method with the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE)Perdew et al. (1996) for the exchange-correlation energy. The valence configurations for Au is [Xe] with cutoff radius 1.32Å; this requires a plane-wave basis set with cutoff energy of 400eV. This gives a lattice constant of 4.171Å for FCC Au, compared with the experimental lattice constant of 4.08Å. Supercells range from 4- through 7-layers of atoms with one Au atom on each layer, and a vacuum gap of 10.5Å to prevent the interaction between slabs under periodic boundary conditions. We use Monkhorst-Pack k-point meshes of for bulk four-atom cells, and for slab supercells; Brillouin-zone integration uses the Methfessel-Paxton methodMethfessel and Paxton (1989) with for electronic occupancies, and the total energy extrapolated to . We represent the charge density and compute energy densities on a grid increasing from to for 4- to 7-layer supercells. Geometry is optimized to reduce forces below 5meV/Å. This gives surface energies of 53.3, 52.7, 52.5 to 52.3meV/Å; this agrees with Zólyomi et al.’s valueZólyomi et al. (2009) of 54meV/Å.

Figure 5 shows the energy change from bulk value for each layer of atoms. Energy density integration of left (or right) 2 layers of 4-layer slab gives a surface energy of meV/Å; 5-, 6-, and 7-layer slabs all give a surface energy of meV/Å. Although the calculated surface energy of 4-layer surface energy is close to the 5-, 6- and 7-layer’s values, the atomic energy distribution shows the center layers of 4-layer slab have not reached the bulk-like behavior. Unlike the convergence test of the traditional total energy calculation, the required thickness of a slab can be directly determined by observing the variation of atomic energy from the bulk value.

### iii.4 HCP Ti with O interstitial

Finally, we consider the formation energy of an oxygen interstitial in the octahedral site of HCP titanium. The formation energy is

 ETi-OI=E(Ti + Oi)−E(Ti)−12E(O2), (31)

where and are the total energy of relaxed supercells with and without an oxygen atom, and is the total energy of oxygen molecule. Our DFT calculations are performed with the PAW method with the GGA-PW91 for the exchange-correlation energy. The valence configurations for Ti is [Ne] with cutoff radius 1.22Å, and O is [He] with cutoff radius 0.80Å; this requires a plane-wave basis set with cutoff energy of 500eV. This gives a lattice constant of 2.933Å, 4.638Å, and 1.581 for HCP Ti, compared with the experimental lattice constant of =2.951Å, 4.684Å, and 1.587.Wood (1962) The supercell contains 96 Ti atoms () and 1 O atom. We use a Monkhorst-Pack k-point mesh; Brillouin-zone integration uses the Methfessel-Paxton method with for electronic occupancies, and the total energy extrapolated to . We represent the charge density and compute energy densities on a grid of . Geometry is optimized to reduce forces below 20meV/Å. This gives an oxygen interstitial formation energy of –6.19eV, with a nearest-neighbor distance between Ti and O of 2.08Å. Hennig et al.Hennig et al. (2005) performed GGA-PW91 calculations using ultrasoft Vanderbilt-typeVanderbilt (1990) pseudopotentials in the same supercell to get a formation energy of –6.12eV and nearest-neighbor distances of 2.06–2.09Å.

Figure 6 shows the calculated Ti atomic energy change from bulk value for each shell. The change in energy shows the differences from bulk are mainly confined to the first two shells, with bulk-like behavior for shells further away from the oxygen atom. The energy density and charge density oscillates and decays away from the interstitial. They peak at the 6th shell and the 13th shell with a wavelength of 1.9Å. Weiss’s Compton profileWeiss (1973) measured the Fermi momentum of Ti as , which corresponds to a Friedel oscillation wavelength of 1.5Å. Jepson’sJepsen (1975) earlier calculation using linear muffin-tin-orbital method obtained the Fermi energy of Ti as 0.667Ryd which corresponds to the Friedel oscillation wavelength of 2.0Å. Adding the atomic energy change of first two shells from the O interstitial and the atomic energy change of O atom, we obtain the interstitial formation energy of eV, which agrees with the total-energy calculation.

## Iv Conclusions

We implement the energy density method for PAW and USPPs for the planewave DFT code vasp; and analyze surface energies from the energy density in the surface region; and vacancy and interstitial formation energies from the energy density in the point defect region. The method can be applied to surfaces and defects in a variety of system, and produces defect formation energies comparable to well-converged total energy calculations. Furthermore, the energy density determines the distribution of energy near the defect or surface, and the sufficiency of a supercell without a separate convergence test. It can also give separate defect formation energies from a single supercell calculation.

###### Acknowledgements.
This research was supported by NSF under grant number DMR-1006077 and through the Materials Computation Center at UIUC, NSF DMR-0325939, and with computational resources from NSF/TeraGrid provided by NCSA and TACC.

## References

• Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. A, 140, 1133 (1965).
• Chetty and Martin (1992) N. Chetty and R. M. Martin, Phys. Rev. B, 45, 6074 (1992a).
• Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B, 43, 1993 (1991).
• Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 45, 566 (1980).
• Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B, 23, 5048 (1981).
• Ewald (1921) P. P. Ewald, Ann. Phys., 369, 253 (1921).
• Chetty and Martin (1992) N. Chetty and R. M. Martin, Phys. Rev. B, 45, 6089 (1992b).
• Zhang and Wei (2004) S. B. Zhang and S. H. Wei, Phys. Rev. Lett., 92, 086102 (2004).
• Rapcewicz et al. (1998) K. Rapcewicz, B. Chen, B. Yakobson,  and J. Bernholc, Phys. Rev. B, 57, 7281 (1998).
• Ramprasad (2002) R. Ramprasad, J. Phys.: Condens. Matter, 14, 5497 (2002).
• Blöchl (1994) P. E. Blöchl, Phys. Rev. B, 50, 17953 (1994).
• Vanderbilt (1990) D. Vanderbilt, Phys. Rev. B, 41, 7892 (1990).
• Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B, 54, 11169 (1996).
• Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B, 59, 1758 (1999).
• Bader (1990) R. F. Bader, Atoms in Molecules: A Quantum Theory (Oxford University Press: Oxford, 1990).
• Rodríguez et al. (2009) J. I. Rodríguez, P. W. Ayers, A. W. Götz,  and F. L. Castillo-Alvarado, J. Chem. Phys., 131, 021101 (2009).
• Martin (2004) R. M. Martin, Electronic structure: Basic Theory and Practical Methods (Cambridge University Press, 2004).
• Slater (1937) J. C. Slater, Phys. Rev., 51, 846 (1937).
• Yu and Trinkle (2010) M. Yu and D. R. Trinkle, arXiv:1010.4916 (2010).
• Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 13, 5188 (1976).
• Moll et al. (1996) N. Moll, A. Kley, E. Pehlke,  and M. Scheffler, Phys. Rev. B, 54, 8844 (1996).
• Qian et al. (1988) G. X. Qian, R. M. Martin,  and D. J. Chadi, Phys. Rev. B, 37, 1303 (1988).
• Choudhury et al. (2008) R. Choudhury, D. R. Bowler,  and M. J. Gillan, J. Phys.: Condens. Matter, 20, 235227 (2008).
• Messmer and Bilello (1981) C. Messmer and J. C. Bilello, J. Appl. Phys., 52, 4623 (1981).
• Wright (2006) A. F. Wright, Phys. Rev. B, 74, 165116 (2006).
• Perdew et al. (1996) J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett., 77, 3865 (1996).
• Puska et al. (1998) M. J. Puska, S. Pöykkö, M. Pesola,  and R. M. Nieminen, Phys. Rev. B, 58, 1318 (1998).
• Dannefaer et al. (1986) S. Dannefaer, P. Mascher,  and D. Kerr, Phys. Rev. Lett., 56, 2195 (1986).
• Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B, 45, 13244 (1992).
• Kane (1985) E. O. Kane, Phys. Rev. B, 31, 5199 (1985).
• McMurry et al. (1967) H. L. McMurry, A. W. Solbrig Jr.,  and J. K. Boyter, J. Phys. Chem. Solids, 28, 2359 (1967).
• Methfessel and Paxton (1989) M. Methfessel and A. T. Paxton, Phys. Rev. B, 40, 3616 (1989).
• Zólyomi et al. (2009) V. Zólyomi, L. Vitos, S. K. Kwon,  and J. Kollár, J. Phys.: Condens. Matter, 21, 095007 (2009).
• Wood (1962) R. M. Wood, Proc. Phys. Soc., 80, 783 (1962).
• Hennig et al. (2005) R. G. Hennig, D. R. Trinkle, J. Bouchet, S. G. Srinivasan, R. C. Albers,  and J. W. Wilkins, Nat. Mater., 4, 129 (2005).
• Weiss (1973) R. J. Weiss, Philos. Mag., 27, 1461 (1973).
• Jepsen (1975) O. Jepsen, Phys. Rev. B, 12, 2988 (1975).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters

316645

How to quickly get a good answer:
• Keep your question short and to the point
• Check for grammar or spelling errors.
• Phrase it like a question
Test
Test description