Total Energy

The total energy of a system is defined as the sum of the kinetic energyE_kin, the electrostatic energyE_es, the exchange-correlation contributionE_xc, and the external energyE_ext

E_tot =E_kin +E_es+E_xc+E_ext. (3.32) The external energy is a consequence pf an applied external fieldV_ext^σ (r)and, hence, its contribution to the total energy is given as

E_ext =X

σ

Z

d³rnσ(r)V_ext^σ (r). (3.33)

The total kinetic energy E_kin can be evaluated directly from the auxiliary Kohn-Sham states

E_kin =X

nσk

f_nσkhΨ_nσk|Tˆ|Ψ_nσki (3.34)

3.4. Total Energy 35 with the occupation numbersf_nσkor indirectly via the eigenvalue sum and a dou-ble counting correction as

E_kin =X

nσk

f_nσkǫ_nσk−X

σ

Z

d³rnσ(r)V_eff^σ(r). (3.35) The exchange-correlation energy in the local-density approximation reads

E_xc = Z

d³rn(r)ǫ_xc[n](r) (3.36)

whereǫ_xc[n]is the exchange-correlation energy density of the homogeneous elec-tron gas at densityn. The electrostatic energy contributionE_esis discussed in detail in the following section.

3.4.1. The Electrostatic Energy

The electrostatic contribution to the total energy consists of three terms,

E_es=E_ZZ+E_eZ+E_ee (3.37)

the core-core interaction, E_ZZ, the electron-core interaction, E_eZ, and the electron-electron term,E_ee, i.e. the electrostatic self-interaction of an electron densityn(r). The latter is given as

E_ee = 1 where the electrostatic potential of a charge distributionρ(r)is given by

V_es[ρ](r) = Z

d³r^′ ρ(r^′)

|r−r^′|. (3.40)

The evaluation of the electrostatic potential of a charge distribution that is not neutral in the spatial average is difficult. It involves either boundary values in the case of isolated systems or it requires the tuning for the convergence of an Ewald summation technique for periodic boundary conditions [71]. To overcome this dif-ficulty, we consider the electrostatics of the atomic cores. The attractive Coulomb potentials of the atomic nuclei at positionsR^a create an external potential which is experienced by the electrons. The contribution to the electrostatic energy is thus given by where an additional Ewald summation over the atoms is required in the case of periodic boundary conditions. The charge distribution associated with the atomic nuclei is equivalent toZ^aδ(r−R^a). Z^ais the number of unit charges in the nucleus i.e. a negative integer. Finally, the core-core interaction is given by

E_ZZ= 1

Note that while each core does not interact electrostatically with itself it does inter-act with its periodic images, if applicable.

We can simplify the difficulties associated to Ewald summations of charged sys-tems by introducing a generalized charge density of electrons and protons. If the

3.4. Total Energy 37 total number of charges (electrons minus protons) vanishes, then the three con-tributions to the electrostatic energy can be written as electrostatic energy of the neutral charge densityρ(r)defined by

ρ(r) =n(r) +X

where the prime onρreminds us of the omission of the self-interaction of the cores.

3.4.2. Electrostatics and the Pseudocharge Construction

In analogy to the PAW transformed expectation values introduced above we split the electrostatic potentialV_es[ρ]defined in Equation (3.40) and the electrostatic en-ergy yielding

The smooth contribution ˜V_es(r)can be represented with a reasonable number of basis functions, e.g. plane waves or grid points, that do not depend on the atomic positions. The two contributions V_es^a, ˜V_es^a are treated on a radial grid inside the augmentation spheresS^aaround each atoma. The radial grids are limited by the radius of the augmentation spheresr^a_aug.

The electrostatic potential of a charge distribution is found by solving the Poisson equation

∆_rV_es[ρ](r) = −4πρ(r) (3.47) with appropriate boundary conditions. In the case of finite systems, we align the electrostatic potential such that it assumes zero for large distances. Due to the treat-ment of neutral charge distributions the potential approaches zero as fast as∼ ^D·^r

|r|³ , whereDis the total dipole moment, or even faster ifDvanishes. In the special case of three periodic boundaries, the value ofV_es(r)may be shifted by an arbitrary con-stant. To fix this, we demand that the spatial average of the electrostatic potential vanishes,R

d³rV_es(r) = 0.

The solution of Poisson’s Equation (3.47) can be found according to Jackson [72]

via the Green function

G_es(r,r^′) = 1

|r−r^′| (3.48)

since

∆_r^′G_es(r,r^′) =∆_r^′ 1

|r−r^′| = −4πδ(r−r^′) (3.49) in Hartree atomic units.

The integration overr^′shows that the electrostatic potential depends onρ(r^′)in all space. For this reason simply replacingρ(r^′) by a smooth quantity inside the spheres would result in a different electrostatic potential in the interstitial region and is thus not suitable. However, we can expand the Coulomb kernel|r−r^′|⁻¹in terms of multipoles such that

1

We then splitρ(r)into three parts, a smooth contribution ˜ρ(r)and two contribu-tionsρ^a(r), ˜ρ^a(r)treated on the radial grid, i.e.

ρ(r) =ρ(r) +˜ X

a

[ρ^a(r) −ρ˜^a(r)]. (3.51) However, in this situation the non-locality of the electrostatic problem requires a special treatment. We demand that ˜ρ(r) reproduces the same electrostatic multi-pole moments asρ(r)inside each sphere. Hence, the atomic correctionρ^a(r)-˜ρ^a(r) must have only vanishing multipole moments. Since ˜ρ(r)coincides with ρ(r) in the interstitial region, we need to find the two identical charge distributions ˜ρ(r) and ˜ρ^a(r)that fulfill this requirement.

Again, ˜ρ(r)needs to be smooth enough to be accurately represented, so it can neither feature the rapid oscillations of the atomic core densities nor the nuclear chargesZ^aδ(|r−R^a|). At this point we introduce the concept of the pseudocharge construction. The compensator function ˆg^a_L(r−R^a)is strictly localized inside the sphereS^a, centered atR^aand has unity multipole moment

Z

d³rgˆL(r)|r|^ℓ′YL^′(rˆ) =δLL^′. (3.52) Using compensator functions, we can identify ˜ρ(r)as

˜

3.4. Total Energy 39 such thatV_es[ρ](˜ r)is the correct electrostatic potential outside of each sphere and shows the correct asymptotic behavior for |r| → ∞. The summation over ℓ is stopped at 2ℓ_max, where ℓ_max is the ℓ-cutoff of projectors and partial waves. The multipole moments∆q^a_L in the implementation are stored ina%qlm(1:).

3.4.3. Electrostatic Potential in the Spheres

For notational convenience, we assumeR^a=0 in this section.

As we can see from Equation (3.50), the Coulomb kernel |r−r^′|⁻¹ is diagonal with respect to L in an ℓm-representation. Given that all densities are stored on radial grid asρ(r,L), the electrostatic problem reduces to the radial integrals

V_es[ρ(r₂,L)](r₁,L) = radiir₂to findV_es(r₁,L). In a real calculation, the second integral is only executed up to the boundary of the sphereS^a.

To find the correct electrostatic potentials and electrostatic energy contributions in the spheres, we start by generating a preliminary potential, ˜Ves^a,pre[ρ˜^a(r^′,L)](r,L) andV_es^a,pre[ρ^a(r^′,L)](r,L), in the sphere from the true densityρ^a and smooth den-sity ˜ρ^a, where we disregard the boundary conditions for both.

Also, due to the non-locality of the electrostatic problem, we cannot findV_es^a(r,L) fromn^a(r,L)and ˜V_es^a(r,L)from ˜n^a(r,L)alone, respectively. Even if ˜n^a(r,L)matches (3.52)), such that we can multiply Equation (3.58) with the compensators ˆg^a_L and

integrate over space. Then where the first term is an integral of the compensators with the smooth electro-static potential on the grid (c.f. a%vlm(1:)) that contains information about exter-nal potential shifts, electrical fields and higher multipoles. In contrast to∆q^a_L (c.f.

a%qlm(1:)) which corrects the multipole moment for the far field, the shifts v^a_L carry information about the electrostatics outside the spheres into the sphere and generate the correct near field

V˜_es^a[ρ˜^a(r^′,L)](r,L) =V˜_es^a,pre[ρ˜^a(r^′,L)](r,L) +v^a_Lr^ℓ, (3.60) V_es^a[ρ^a(r^′,L)](r,L) =V_es^a,pre[ρ^a(r^′,L)](r,L) +v^a_Lr^ℓ. (3.61)

3.4.4. Electrostatic Energy in the Spheres

We can now compute the electrostatic energy in the spheres. At this point, the treatment of the smooth quantities and the true quantities differs slightly. The electrostatic energy of the smooth density in the sphere is a pure Hartree energy

E˜^a_es = 1 In return, the true electrostatic energyE^a_esin the sphere consists of a Hartree energy and the Coulomb energy of the electron density in the singular potential of the atomic nucleus. It cannot be evaluated byR

V_es^a[ρ^a]ρ^awithρ^a=n^a+Z^aδbecause that would involve a self-interaction of the nucleus. Instead

E^a_es = 1

This way, the electrostatic interaction of a nucleus with itself is excluded.

3.4.5. Compensator Functions

The pseudocharge compensators need to be rather smooth and strictly localized.

We therefore adopt the polynomial construction as described by Weinertet al.[73]

as depicted in Figure 3.5.

ˆ

g^[N_ℓ ^W](r)∼ |r|^ℓ 1− |r|² R²_aug

!(N_W−ℓ)

Yℓm(ˆr) for |r|< R_aug. (3.64)

3.4. Total Energy 41

Figure 3.5.:Strictly localized compensator functions of the Weinert type.

The exponentN_Wcontrols the smoothness of the functions. The main graph shows theℓ=0-compensator for different exponents, whereas the inlet shows the normalized radial functions for various values ofℓ.

3.4.6. Exchange-Correlation

A major advantage of the real-space representation is the ability to directly eval-uate the exchange-correlation energy densityǫ_xc[n](r)and the potentialV_xc[n](r) in a local parametrization as LDA or GGA. However, the true densityndiffers in-side the spheres from the smooth density ˜n. The smooth density ˜n(r)is stored on the real-space grid, and the true densities inside the spheres are stored in(r,L) -representation. Therefore, the exchange-correlation energy has three contributions

E_xc =

Since those exchange-correlation parametrizations known to work are non-linear, the evaluation ofǫ_xc and the matrix elements of partial waves withV_xc has to be performed fully numerically at each step. The densities inside the spheres are

eas-ily represented in a radial representation and spherical harmonics for the angular part, i.e. ˜n^a(r,L), where r will be sampled from the origin to the sphere radius in a usually non-linear fashion. This allows us to capture the rapid oscillations in the all-electron densities with a reasonable number of radial support points. The advantages of an ℓm-representation (or L) are the simplicity of the cutoff param-eter ℓ_max that controls the quality of this representation and the straightforward evaluation of electrostatics. However, the evaluation of the non-linear exchange-correlation potentialǫ_xcandV_xccannot be executed in theℓm-basis. Here, we need to transform into representation that samples the full solid angleΩ. This transfor-mation is performed on a Legendre grid, where we sampleϕandϑin an equidis-tant fashion.

The number of points for this sampling is usually 2M²whereMis chosen>ℓ_max. The transformation then reads

n^a(r,ϑj,ϕi) = X

L

n^a(r,L)YL(ϑj,ϕi), (3.69) and the reverse transformation applied to the potentialV_xc(and alsoǫ_xc) reads

V_xc^a(r,L) =

2πZ

0

dϕ Zπ

0

dϑsinϑ V_xc^a(r,ϑ,ϕ)YL(ϑ,ϕ) (3.70)

≈ XM

j

wj 2M

X

i

V_xc^a(r,ϑj,ϕi)YL(ϑ_j,ϕi)

where the indicesiandjsampleϕandϑ, respectively. A Legendre grid forM=12 is depicted in Figure 3.6. The closer the points are to one of the poles (ϑ=0 orϑ=π), the denser is this grid. Therefore, the weightswjwith

wj =sin(ϑj) (3.71)

need to be included in the reverse transformation.

Im Dokument Real-Space Finite-Difference (Seite 42-51)