The Projector Augmented-wave Method

Carsten Rostgaard
February 3, 2010

Abstract
The purpose of this text is to give a self-contained description of the
basic theory of the projector augmented-wave (PAW) method, as well as
most of the details required to make the method work in practice. These
two topics are covered in the first two sections, while the last is dedicated
to examples of how to apply the PAW transformation when extracting
non-standard quantities from a density-functional theory (DFT) calcula-
tion.
The formalism is based on Blöchl’s original formulation of PAW [1],
and the notation and extensions follow those used and implemented in
the gpaw[2] code.

Contents
1 Formalism 3
1.1 The Transformation Operator . . . . . . . . . . . . . . . . . . . . 3
1.2 The Frozen Core Approximation . . . . . . . . . . . . . . . . . . 6
1.3 Expectation Values . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Total Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.1 The Semi-local Contributions . . . . . . . . . . . . . . . . 9
1.5.2 The Nonlocal Contributions . . . . . . . . . . . . . . . . . 10
1.5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6 The Transformed Kohn-Sham Equation . . . . . . . . . . . . . . 14
1.6.1 Orthogonality . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.6.2 The Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . 15
1.7 Forces in PAW . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Implementing PAW 17
2.1 Atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.1 The Radial Kohn-Sham Equation . . . . . . . . . . . . . . 19
2.2 The Atomic Data Set of PAW . . . . . . . . . . . . . . . . . . . . 19

3 Non-standard Quantities 23
3.1 The External Potential in PAW . . . . . . . . . . . . . . . . . . . 23
3.2 All-electron Density . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Wannier Orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 Local Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.1 Projected Density of States . . . . . . . . . . . . . . . . . 25
3.4.2 Local Magnetic Moments . . . . . . . . . . . . . . . . . . 26
3.4.3 LDA + U . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5 Coulomb Integrals . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.1 Exact Exchange . . . . . . . . . . . . . . . . . . . . . . . 28
3.5.2 Optimized Effective Potential . . . . . . . . . . . . . . . . 30

1
References 31

2
1 Formalism
By the requirement of orthogonality, DFT wave functions have very sharp fea-
tures close to the nuclei, as all the states are non-zero in this region. Further
out only the valence states are non-zero, resulting in much smoother wave func-
tions in this interstitial region. The oscillatory behavior in the core regions,
requires a very large set of plane waves, or equivalently a very fine grid, to be
described correctly. One way of solving this problem is the use of pseudopoten-
tials in which the collective system of nuclei and core electrons are described
by an effective, much smoother, potential. The KS equations are then solved
for the valence electrons only. The pseudopotentials are constructed such that
the correct scattering potential is obtained beyond a certain radius from the
core. This method reduces the number of wave functions to be calculated, since
the pseudo potentials only have to be calculated and tabulated once for each
type of atom, so that only calculations on the valence states are needed. It
justifies the neglect of relativistic effects in the KS equations, since the valence
electrons are non-relativistic (the pseudopotentials describing core states are of
course constructed with full consideration of relativistic effects). The technique
also removes the unwanted singular behavior of the ionic potential at the lattice
points.
The drawback of the method is that all information on the full wave function
close to the nuclei is lost. This can influence the calculation of certain properties,
such as hyperfine parameters, and electric field gradients. Another problem is
that one has no before hand knowledge of when the approximation yields reliable
results.
A different approach is the augmented-plane-wave method (APW), in which
space is divided into atom-centered augmentation spheres inside which the wave
functions are taken as some atom-like partial waves, and a bonding region out-
side the spheres, where some envelope functions are defined. The partial waves
and envelope functions are then matched at the boundaries of the spheres.
A more general approach is the projector augmented wave method (PAW)
presented here, which offers APW as a special case[3], and the pseudopotential
method as a well defined approximation[4]. The PAW method was first proposed
by Blöchl in 1994[1].

1.1 The Transformation Operator

The features of the wave functions, are very different in different regions of
space. In the bonding region it is smooth, but near the nuclei it displays rapid
oscillations, which are very demanding on the numerical representation of the
wave functions. To address this problem, we seek a linear transformation T̂
which takes us from an auxiliary smooth wave function |ψ̃n ⟩ to the true all
electron Kohn-Sham single particle wave function |ψn ⟩

|ψn ⟩ = T̂ |ψ̃n ⟩, (1)

where n is the quantum state label, containing a k index, a band index, and a
spin index.
This transformation yields the transformed KS equations

T̂ † H
b T̂ |ψ̃n ⟩ = ϵn T̂ † T̂ |ψ̃n ⟩, (2)

3
which needs to be solved instead of the usual KS equation. Now we need to
define T̂ in a suitable way, such that the auxiliary wave functions obtained from
solving 2 becomes smooth.
Since the true wave functions are already smooth at a certain minimum
distance from the core, T̂ should only modify the wave function close to the
nuclei. We thus define X
T̂ = 1 + T̂ a , (3)
a

where a is an atom index, and the atom-centered transformation, T̂ a , has no

effect outside a certain atom-specific augmentation region |r − Ra | < rca . The
cut-off radii, rca should be chosen such that there is no overlap of the augmen-
tation spheres.
Inside the augmentation spheres, we expand the true wave function in the
partial waves ϕai , and for each of these partial waves, we define a corresponding
auxiliary smooth partial wave ϕ̃ai , and require that

|ϕai ⟩ = (1 + T̂ a )|ϕ̃ai ⟩ ⇔ T̂ a |ϕ̃ai ⟩ = |ϕai ⟩ − |ϕ̃ai ⟩ (4)

for all i, a. This completely defines T̂ , given ϕ and ϕ̃.

Since T̂ a should do nothing outside the augmentation sphere, we see from 4
that we must require the partial wave and its smooth counterpart to be identical
outside the augmentation sphere

∀a, ϕai (r) = ϕ̃ai (r), for r > rca

where ϕai (r) = ⟨r|ϕai ⟩ and likewise for ϕ̃ai .

If the smooth partial waves form a complete set inside the augmentation
sphere, we can formally expand the smooth all electron wave functions as
X
a
|ψ̃n ⟩ = Pni |ϕ̃ai ⟩, for |r − Ra | < rca (5)
i

a
where Pni are some, for now, undetermined expansion coefficients.
Since |ϕai ⟩ = T̂ |ϕ̃ai ⟩ we see that the expansion
X
a
|ψn ⟩ = T̂ |ψ̃n ⟩ = Pni |ϕai ⟩, for |r − Ra | < rca (6)
i

a
has identical expansion coefficients, Pni .
a
As we require T̂ to be linear, the coefficients Pni must be linear functionals
of |ψ̃n ⟩, i.e. Z
a
Pni = ⟨p̃ai |ψ̃n ⟩ = drp̃a∗ a
i (r − R )ψ̃n (r), (7)

where |p̃ai ⟩ are some fixed functions termed smooth projector functions.
As there is no overlap between the augmentation spheres, we P expect the one
center expansion of the smooth all electron wave function, |ψ̃na ⟩ = i |ϕ̃ai ⟩⟨p̃ai |ψ̃n ⟩
to reduce to |ψ̃n ⟩ itself inside the augmentation sphere defined by a. Thus, the
smooth projector functions must satisfy
X
|ϕ̃ai ⟩⟨p̃ai | = 1 (8)
i

4
inside each augmentation sphere. This also implies that

⟨p̃ai1 |ϕ̃ai2 ⟩ = δi1 ,i2 , for |r − Ra | < rca (9)

i.e. the projector functions should be orthonormal to the smooth partial waves
inside the augmentation sphere. There are no restrictions on p̃ai outside the
augmentation spheres, so for convenience we might as well define them as local
functions, i.e. p̃ai (r) = 0 for r > rca .
Note that the completeness relation 8 is equivalent to the requirement that
p̃ai should produce the correct expansion coefficients of 5-6, while 9 is merely
an implication of this restriction. Translating 8 to an explicit restriction on the
projector functions is a rather involved procedure, but according to Blöchl, [1],
the most general form of the projector functions is:
X
⟨p̃ai | = ({⟨fka |ϕ̃al ⟩})−1 a
ij ⟨fj |, (10)
j

where |fja ⟩ is any set of linearly independent functions. The projector functions
will be localized if the functions |fja ⟩ are localized.
Using the completeness relation 8, we see that
X X
T̂ a = T̂ a |ϕ̃ai ⟩⟨p̃ai | = |ϕai ⟩ − |ϕ̃ai ⟩ ⟨p̃ai |,


i i

where the first equality is true in all of space, since 8 holds inside the augmen-
tation spheres and outside T̂ a is zero, so anything can be multiplied with it.
The second equality is due to 4 (remember that |ϕai ⟩ − |ϕ̃ai ⟩ = 0 outside the
augmentation sphere). Thus we conclude that
XX
|ϕai ⟩ − |ϕ̃ai ⟩ ⟨p̃ai |.


T̂ = 1 + (11)
a i

To summarize, we obtain the all electron KS wave function ψn (r) = ⟨r|ψn ⟩

from the transformation
XX
ϕai (r) − ϕ̃ai (r) ⟨p̃ai |ψ̃n ⟩,


ψn (r) = ψ̃n (r) + (12)
a i

where the smooth (and thereby numerically convenient) auxiliary wave function
ψ̃n (r) is obtained by solving the eigenvalue equation 2.
The transformation 12 is expressed in terms of the three components: a)
the partial waves ϕai (r), b) the smooth partial waves ϕ̃ai (r), and c) the smooth
projector functions p̃ai (r).
The restriction on the choice of these sets of functions are: a) Since the
partial- and smooth partial wave functions are used to expand the all electron
wave functions, i.e. are used as atom-specific basis sets, these must be complete
(inside the augmentation spheres). b) the smooth projector functions must
satisfy 8, i.e. be constructed according to 10. All remaining degrees of freedom
are used to make the expansions converge as fast as possible, and to make the
functions termed ‘smooth’, as smooth as possible. For a specific choice of these
sets of functions, see section 2.2. As the partial- and smooth partial waves are
merely used as basis sets they can be chosen as real functions (any imaginary

5
parts of the functions they expand, are then introduced through the complex
a
expansion coefficients Pni ). In the remainder of this document ϕ and ϕ̃ will be
assumed real.
Note that the sets of functions needed to define the transformation are sys-
tem independent, and as such they can conveniently be pre-calculated and tab-
ulated for each element of the periodic table.
For future convenience, we also define the one center expansions
X
ψna (r) = ϕai (r)⟨p̃ai |ψ̃n ⟩, (13a)
i
X
ψ̃na (r) = ϕ̃ai (r)⟨p̃ai |ψ̃n ⟩. (13b)
i

In terms of these, the all electron KS wave function is

X
ψna (r − Ra ) − ψ̃na (r − Ra ) .


ψn (r) = ψ̃n (r) + (14)
a

So what have we achieved by this transformation? The trouble of the original

KS wave functions, was that they displayed rapid oscillations in some parts of
space, and smooth behavior in other parts of space. By the decomposition 12 we
have separated the original wave functions into auxiliary wave functions which
are smooth everywhere and a contribution which contains rapid oscillations, but
only contributes in certain, small, areas of space. This decomposition is shown
on the front page for the hydrogen molecule. Having separated the different
types of waves, these can be treated individually. The localized atom centered
parts, are indicated by a superscript a, and can efficiently be represented on
atom centered radial grids. Smooth functions are indicated by a tilde ˜. The
delocalized parts (no superscript a) are all smooth, and can thus be represented
on coarse Fourier- or real space grids.

1.2 The Frozen Core Approximation

In the frozen core approximation, it is assumed that the core states are natu-
rally localized within the augmentation spheres, and that the core states of the
isolated atoms are not changed by the formation of molecules or solids. Thus
the core KS states are identical to the atomic core states

|ψnc ⟩ = |ϕa,core
α ⟩,

where the index n on the left hand site refers to both a specific atom, a, and an
atomic state, α.
In this approximation only valence states are included in the expansions of
|ψn ⟩, 6, and |ψ̃n ⟩, 5.
Figure 1, shows the atomic states of Platinum in its ground state, obtained
with an atomic DFT program at an LDA level, using spherical averaging, i.e. a
spin-compensated calculation, assuming the degenerate occupation 9/10 of all
5d states, and both of the 6s states half filled. It is seen that at the typical length
of atomic interaction (the indicated cut-off rc = 2.5 Bohr is approximately half
the inter-atomic distance in bulk Pt), only the 5d and 6s states are non-zero.
Note that the frozen core approximation is not a prerequisite for PAW. See
e.g. [5] for a relaxation of this requirement.

6
 Figure 1: The core states of Platinum.

1.3 Expectation Values

The expectation value of an operator O
b is, within the frozen core approximation,
given by
val
X XX core
⟨O⟩
b = fn ⟨ψn |O|ψ
b n⟩ + ⟨ϕa,core
α
b a,core
|O|ϕ α ⟩. (15)
n a α

b n ⟩ = ⟨ψ̃n |T̂ † O
Using that ⟨ψn |O|ψ b T̂ |ψ̃n ⟩, and skipping the state index for nota-
tional convenience, we see that
X X
⟨ψ|O|ψ⟩
b = ⟨ψ̃ + (ψ a − ψ̃ a )|O|
b ψ̃ + (ψ a − ψ̃ a )⟩
a a
X 
b a′ − ψ̃ a′ ⟩ +
X
= ⟨ψ̃|O|
b ψ̃⟩ + ⟨ψ a − ψ̃ a |O|ψ b a − ψ̃ a ⟩ + ⟨ψ a − ψ̃ a |O|
⟨ψ̃|O|ψ b ψ̃⟩
aa′ a
X 
= ⟨ψ̃|O|
b ψ̃⟩ + ⟨ψ a |O|ψ
b a ⟩ − ⟨ψ̃ a |O|
b ψ̃ a ⟩
a
X 
+ ⟨ψ a − ψ̃ a |O|
b ψ̃ − ψ̃ a ⟩ + ⟨ψ̃ − ψ̃ a |O|ψ
b a − ψ̃ a ⟩
a
X ′ ′
+ ⟨ψ a − ψ̃ a |O|ψ
b a − ψ̃ a ⟩.
a̸=a′
(16)

For local operators1 the last two lines does not contribute. The first line, because
|ψ a − ψ̃ a ⟩ is only non-zero inside the spheres, while |ψ̃ − ψ̃ a ⟩ is only non-zero
outside the spheres. The second line simply because |ψ a − ψ̃ a ⟩ is zero outside the
1 Local operator O: b An operator which does not correlate separate parts of space, i.e.
b ′ ⟩ = 0 if r ̸= r′ .
⟨r|O|r

7
spheres, so two such states centered on different nuclei have no overlap (provided
that the augmentation spheres do not overlap).
Reintroducing the partial waves in the one-center expansions, we see that
val
X val
X X X val
X
fn ⟨ψna |O|ψ
b na ⟩ = fn ⟨ϕai1 Pni
a
1
b ai Pni
|O|ϕ 2
a
2
⟩= ⟨ϕai1 |O|ϕ
b ai ⟩
2
a∗ a
fn Pni P ,
1 ni2
n n i1 i2 i1 i2 n
(17)
and likewise for the smooth waves.
Introducing the Hermitian one-center density matrix
X X
Dia1 i2 = a∗ a
fn Pni P =
1 ni2
fn ⟨ψ̃n |p̃ai1 ⟩⟨p̃ai2 |ψ̃n ⟩. (18)
n n

We conclude that for any local operator O,

b the expectation value is

val
X XX  core
XX
⟨O⟩
b = fn ⟨ψ̃n |O|
b ψ̃n ⟩+ ⟨ϕai1 |O|ϕ
b ai ⟩ − ⟨ϕ̃ai |O|
2 1
b ϕ̃ai ⟩ Dia i +
2 1 2
⟨ϕa,core
α
b a,core
|O|ϕ α ⟩.
n a i1 i2 a α
(19)

1.4 Densities
The electron density is obviously a very important quantity in DFT, as all
observables in principle are calculated as functionals of the density. In reality
the kinetic energy is calculated as a functional of the orbitals, and some specific
exchange-correlation functionals also rely on KS-orbitals rather then the density
for their evaluation, but these are still implicit functionals of the density.
To obtain the electron density we need to determine the expectation value
of the real-space projection operator |r⟩⟨r|
X X
n(r) = fn ⟨ψn |r⟩⟨r|ψn ⟩ = fn |ψn (r)|2 , (20)
n n

where fn are the occupation numbers.

As the real-space projection operator is obviously a local operator, we can
use the results 19 of the previous section, and immediately arrive at
val
X XX  core
XX
n(r) = fn |ψ̃n |2 + ϕai1 ϕai2 − ϕ̃ai1 ϕ̃ai2 Dia1 i2 + |ϕa,core
α |2 . (21)
n a i1 i2 a α

To ensure that 21 reproduce the correct density even though some of the
core states are not strictly localized within the augmentation spheres, a smooth
core density, ñc (r), is usually constructed, which is identical to the core density
outside the augmentation sphere, and a smooth continuation inside. Thus the
density is typically evaluated as
X
n(r) = ñ(r) + (na (r) − ña (r)) , (22)
a

8
where
val
X
ñ(r) = fn |ψ̃n (r)|2 + ñc (r) (23a)
n
X
na (r) = Dia1 i2 ϕai1 (r)ϕai2 (r) + nac (r) (23b)
i1 i2
X
ña (r) = Dia1 i2 ϕ̃ai1 (r)ϕ̃ai2 (r) + ñac (r) (23c)
i1 i2

1.5 Total Energies

The total energy of the electronic system is given by:

E[n] = Ts [n] + UH [n] + Vext [n] + Exc [n]. (24)

In this section, the usual energy expression above, is sought re-expressed in terms
of the PAW quantities: the smooth waves and the auxiliary partial waves.
For the local and semi-local functionals, we can utilize 19, while the nonlocal
parts needs more careful consideration.

1.5.1 The Semi-local Contributions

The kinetic energy functional Ts = n fn ⟨ψn | − 21 ∇2 |ψn ⟩ is obviously a (semi-)
P
local functional, so we can apply 19 and immediately arrive at:
X
Ts [{ψn }] = fn ⟨ψn | − 12 ∇2 |ψn ⟩
n
val
X X (25)
1 2
Tca ∆Tia1 i2 Dia1 i2


= fn ⟨ψ̃n | − 2 ∇ |ψ̃n ⟩ + + ,
n a

where
core
X
Tca = ⟨ϕa,core
α |− 12 ∇2 |ϕa,core
α ⟩ and ∆Tia1 i2 = ⟨ϕai1 |− 12 ∇2 |ϕai2 ⟩−⟨ϕ̃ai1 |− 12 ∇2 |ϕ̃ai2 ⟩.
α
(26)
For LDA and GGA type exchange-correlation functionals, Exc is likewise, per
definition, a semi-local functional, such that it can be expressed as
X
Exc [n] = Exc [ñ] + (Exc [na ] − Exc [ña ]) . (27)
a

By virtue of 23b-23c we can consider the atomic corrections as functionals of

the density matrix defined in 18, i.e.
X
a
Exc [n] = Exc [ñ] + ∆Exc [{Dia1 i2 }], (28)
a

where
a
∆Exc [{Dia1 i2 }] = Exc [na ] − Exc [ña ]. (29)

9
1.5.2 The Nonlocal Contributions
The Hartree term is both nonlinear and nonlocal, so more care needs to be taken
when introducing the PAW transformation for this expression.
In the following we will assume that there is no ‘true’ external field, such that
Vext [n] is only due to the static nuclei, i.e. it is a sum of the classical interaction
of the electron density with the static ionic potential, and the electrostatic
energy of the nuclei.
We define the total classical electrostatic energy functional as
1 1 ′
((n)) + (n| a Z a ) + a
|Z a ),
P P
EC [n] = UH [n] + Vext [n] = 2 a̸=a′ (Z
(30)
2
where the notation (f—g) indicates the Coulomb integral

f ∗ (r)g(r′ )
ZZ
(f |g) = drdr′ (31)
|r − r′ |

and I have introduced the short hand notation ((f )) = (f |f ). In 30, Z a (r) is
the charge density of the nucleus at atomic site a, which in the classical point
approximation is given by

Z a (r) = −Z a δ(r − Ra ) (32)

with Z a being the atomic number of the nuclei. As the Hartree energy of a
density with non-zero total charge is numerically inconvenient, we introduce
the charge neutral total density
X
ρ(r) = n(r) + Z a (r) (= nelectrons + nnuclei ). (33)
a

In terms of this, the coulombic energy of the system can be expressed by

′ 1
((n + a Z a ))′
P
EC [n] = UH [ρ] = (34)
2
where the prime indicates that one should remember the self-interaction error of
the nuclei introduced in the Hartree energy of the total density. This correction
is obviously ill defined, and different schemes exist for making this correction.
As it turns out, this correction is handled very naturallyP in the PAW formalism.
For now, we will focus on the term ((ρ)) = ((n + a Z a )). If one where to
directly include the expansion of n(r) according to 22, one would get:

((n + a Z a )) = ((ñ + a na − ña + Z a ))

P P
X ′ ′ ′ X
= ((ñ)) + (na − ña + Z a |na − ña + Z a ) + 2 (ñ|na − ña + Z a ),
aa′ a

where in the last expression, the first term is the Hartree energy of the smooth
electron density, which is numerically problematic because of the nonzero total
charge. The second term contains a double summation over all nuclei, which
would scale badly with system size, and the last term involves integrations
of densities represented on incompatible grids (remember that the one-center
densities are represented on radial grids to capture the oscillatory behavior near

10
the nuclei)2 . This is clearly not a feasible procedure. To correct these problem
we add and subtract some atom centered compensation charges Z̃ a :
P h i ′ ′ ′
((n+ a Z̃ a + a Z a − Z̃ a )) = ((ñ+ a Z̃ a ))+ aa′ (na −ña +Z a −Z̃ a |na −ña +Z a −Z̃ a )
P P P

′
X
(ñ + a′ Z̃ a |na − ña + Z a − Z̃ a ).
P
+2
a

If we define Z̃ (r) in such a way that na (r) − ña (r) + Z a (r) − Z̃ a (r) has no
a

multipole moments, i.e.

Z
drrl YL (r −
d Ra )(na − ña + Z a − Z̃ a ) = 0 (35)

for all a, the potentials of these densities are zero outside their respective aug-
mentation spheres (L = (l, m) is a collective angular- and magnetic quantum
number). Exploiting this feature, the Coulomb integral reduce to
((n + a Z a )) = ((ñ + a Z̃ a )) + a ((na − ña + Z a − Z̃ a )) + 2 a (ña + Z̃ a |na − ña + Z a − Z̃ a )
P P P P
P  
= ((ñ + a Z̃ a )) + a ((na + Z a )) − ((ña + Z̃ a ))
P

where it has been used that inside the augmentation spheres ñ = ña . In this
expression, we have circumvented all of the previous problems. None of the
terms correlates functions on different grids, there is only a single summation
over the atomic sites, and furthermore the onlyPthing that has to be evaluated in
the full space is the Hartree energy of ñ(r) + a Z̃ a (r) which is charge neutral
(see eq. 42).
Inserting the final expression in 30, we see that
1 1 X a 
EC [n] = ((ñ + a Z̃ a )) + ((n + Z a ))′ − ((ña + Z̃ a ))
P
2 2 a
(36)
1 X  
= UH [ρ̃] + ((na )) + 2(na |Z a ) − ((ña + Z̃ a ))
2 a

where we have introduced the smooth total density

X
ρ̃(r) = ñ + Z̃ a (r). (37)
a

Note that the problem with the self interaction error of the nuclei could easily
be resolved once it was moved to the atom centered part, as handling charged
densities is not a problem on radial grids.
To obtain an explicit expression for the compensation charges, we make a
multipole expansion of Z̃ a (r)
X
Z̃ a = QaL g̃L
a
(r), (38)
L
a
where g̃L (r) is any smooth function localized within |r − Ra | < rca , satisfying
Z
drrl YL (r −
d a
Ra )g̃L ′ (r) = δLL′ . (39)

2 One could separate the terms in other ways, but it is impossible to separate the smooth

and the localized terms completely.

11
 Plugging the expansion 38 into equations 35, we see that the expansion
coefficients QaL from must be chosen according to
Z X
QaL = drrl YL (r̂) (na (r) − ña (r) + Z a (r)) = ∆a δl,0 + ∆aLi1 i2 Dia1 i2 (40)
i1 i2

where
a
Z √
∆ = dr (nac (r) − ñac (r)) − Z a / 4π (41a)
Z
∆aLi1 i2 = drrl YL (r̂)[ϕai1 (r)ϕai2 (r) − ϕ̃ai1 (r)ϕ̃ai2 (r)] (41b)

and it has been used that the core densities are spherical nac (r) = nac (r)Y00 (r̂)
(we consider only closed shell frozen cores). This completely defines the com-
pensation charges Z̃ a (r).
Note that the special case l = 0 of 35, implies that
Z h i
dr na − ña + Z a − Z̃ a = 0

⇓
Z " # Z " #
X X
a a
dr ñ(r) + Z̃ (r) = dr n(r) + Z (r)
a a
⇕
Z Z
dr ρ̃(r) = dr ρ(r) = 0 (42)

i.e. that the smooth total density has zero total charge, making the evaluation
of the Hartree energy numerically convenient.
In summary, we conclude that the classical coulomb interaction energy which
′
in the KS formalism was given by EC [n] = UH [ρ], in the PAW formalism be-
comes a pure Hartree energy (no self-interaction correction) of the smooth total
density ρ̃ plus some one-center corrections
X
a
EC [n] = UH [ρ̃] + ∆EC [{Dia1 i2 }] (43)
a

where the corrections

a
∆EC [{Dia1 i2 }] = 21 ((na )) + (na |Z a ) − 12 ((ña )) − (ña |Z̃ a ) − 12 ((Z̃ a ))
na (r) X a a a
Z
1 a a
= 2 [((nc )) − ((ñc ))] − Z a
dr c − QL (ñc |g̃L )
r
L
" #
ϕai1 (r)ϕai2 (r) X a a a a
X Z
a∗ a a a a a a a
+ Di1 i2 (ϕi1 ϕi2 |nc ) − (ϕ̃i1 ϕ̃i2 |ñc ) − Z dr − QL (ϕ̃i1 ϕ̃i2 |g̃L )
i1 i2
r
L
1 X h i 1X a a a a
+ Dia∗
1 i2
(ϕ a a
ϕ |ϕ a a
i1 i2 i3 i4ϕ ) − (ϕ̃ a a
ϕ̃ |ϕ̃ a a
ϕ̃
i1 i2 i3 i4 ) Dia3 i4 − Q Q ′ (g̃ |g̃ ′ )
2i i i i 2 ′ L L L L
1 2 3 4 LL

Using that the potential of a spherical harmonic (times some radial function) is
a a
itself a spheical harmonic of the same angular momentum, we see that (g̃L |g̃L′ ) ∝

12
δLL′ and (ñac |g̃L
a
) ∝ δL0 . Noting that QaL by virtue of 40 is a functional of the
density matrix, and inserting this, we get
X X
a
∆EC = ∆C a + ∆Cia1 i2 Dia1 i2 + Dia∗
1 i2
∆Cia1 i2 i3 i4 Dia3 i4 (44)
i1 i2 i1 i2 i3 i4

where

a
√ Z
nac (r)
= 12 ((nac )) ((ñac )) a 2 a a
(ñac |g̃00
a a
 
∆C − − (∆ ) ((g̃00 )) −∆ ) − 4πZ dr
r
(45)
ϕai1 (r)ϕai2 (r)
Z
∆Cia1 i2 =(ϕai1 ϕai2 |nac ) − (ϕ̃ai1 ϕ̃ai2 |ñac ) − Z a dr
r
 
a a a a a a a a a
− ∆ (ϕ̃i1 ϕ̃i2 |g̃00 ) − ∆00,i1 i2 ∆ (ñc |g̃00 ) + ((g̃00 )) (46)
h i
∆Cia1 i2 i3 i4 = 21 (ϕai1 ϕai2 |ϕai3 ϕai4 ) − (ϕ̃ai1 ϕ̃ai2 |ϕ̃ai3 ϕ̃ai4 )
Xh i
− 12 ∆aLi1 i2 (ϕ̃ai1 ϕ̃ai2 |g̃L
a
) + ∆aLi3 i4 (ϕ̃ai3 ϕ̃ai4 |g̃L
a
) + ∆aLi1 i2 ((g̃L
a
))∆aLi3 i4
L
(47)

Note that all integrals can be limited to the inside of the augmentation sphere.
For example (ϕai1 ϕai2 |nac ) has contributions outside the augmentation sphere,
but these are exactly canceled by the contributions outside the spheres of
(ϕ̃ai1 ϕ̃ai2 |ñac ), in which region the two expressions are identical.
The ∆Cia1 i2 i3 i4 tensor has been written in a symmetric form, such that it is
invariant under the following symmetry operations:

i1 ↔ i2 i3 ↔ i4 i1 i2 ↔ i3 i4 (48)

1.5.3 Summary
Summing up all the energy contributions, we see that the Kohn-Sham total
energy
′
E[n] = Ts [{ψn }] + UH [ρ] + Exc [n]
can be separated into a part calculated on smooth functions, Ẽ, and some atomic
corrections, ∆E a , involving quantities localized around the nuclei only.
X
E = Ẽ + ∆E a (49)
a

where the smooth part

Ẽ = Ts [{ψ̃n }] + UH [ρ̃] + Exc [ñ] (50)

is the usual energy functional, but evaluated on the smooth functions ñ and ρ̃
instead of n and ρ, and with the soft compensation charges Z̃ a instead of the
nuclei charges Z a (r). The corrections are given by
X
X
∆E a = ∆Tca +∆C a + ∆Tia1 i2 + ∆Cia1 i2 + Dia∗
1 i2
∆Cia1 i2 i3 i4 Dia3 i4 +∆Exc
a
[{Dia1 i2 }]
i1 i2 i1 i2 i3 i4
(51)

13
where Tca , ∆Tia1 i2 , ∆Cia1 i2 , and ∆Cia1 i2 i3 i4 are system independent tensors that
can be pre-calculated and stored for each specie in the periodic table of elements.
Both the Hamiltonian and the forces can be derived from the total energy
functional 49 as will be shown in the following two sections.

1.6 The Transformed Kohn-Sham Equation

The variational quantity in the PAW formalism is the smooth wave function ψ̃n .
From this, all other quantities, such as the density matrix, the soft compensation
charges, the transformation operator, etc. are determined by various projections
of ψ̃n onto the projector functions, and expansions in our chosen basis functions,
the partial and smooth partial waves. To obtain the smooth wave functions, we
need to solve the eigenvalue equation

H
e ψ̃n (r) = ϵn Ŝ ψ̃n (r), (52)
b

where the overlap operator Ŝ = T̂ † T̂ and H = T̂ † H

b T̂ is the transformed
be
Hamiltonian.

1.6.1 Orthogonality
In the original form, the eigen states of the KS equation where orthogonal, i.e.
⟨ψn |ψm ⟩ = δnm while in the transformed version

⟨ψ̃n |T̂ † T̂ |ψ̃m ⟩ = δnm (53)

i.e. the smooth wave function are only orthogonal with respect to the weight Ŝ.
The explicit form of the overlap operator is

Ŝ = T̂ † T̂
 †  
= 1 + a T̂ a 1 + a T̂ a
P P
X 
=1+ T̂ a† + T̂ a + T̂ a† T̂ a
a
XhX X X X
=1+ |p̃ai1 ⟩(⟨ϕai1 | − ⟨ϕ̃ai1 |) |ϕ̃ai2 ⟩⟨p̃ai2 | + |ϕ̃ai2 ⟩⟨p̃ai2 | (|ϕai1 ⟩ − |ϕ̃ai1 ⟩)⟨p̃ai1 |
a i1 i2 i2 i1
X X i
+ |p̃ai1 ⟩(⟨ϕai1 | − ⟨ϕ̃ai1 |) (|ϕai2 ⟩ − |ϕ̃ai2 ⟩)⟨p̃ai2 |
i1 i2
XX
=1+ |p̃ai1 ⟩(⟨ϕai1 |ϕai2 ⟩ − ⟨ϕ̃ai1 |ϕ̃ai2 ⟩)⟨p̃ai2 |
a i1 i2
XX √
=1+ |p̃ai1 ⟩ 4π∆a00,i1 i2 ⟨p̃ai2 |
a i1 i2
(54)

The orthogonality condition 53 must be kept in mind when applying numerical

schemes for solving 52. For example plane waves are no longer orthogonal, in
the sense that ⟨G|Ŝ|G′ ⟩ =
̸ δG,G′ .

14
1.6.2 The Hamiltonian
To determine the transformed Hamiltonian, one could apply the transformation
H = T̂ † H
b T̂ directly, which would be straight forward for the local parts of H,
be b
but to take advantage of the trick used to determine the total energy of the
nonlocal term (EC [n]), we make use of the relation

δE
= fn H ψ̃n (r). (55)
be
∗
δ ψ̃n (r)

Using this, we get

δE δ h i
= Ts [{ψ̃n }] + Exc [ñ] + UH [ρ̃] + ∆E a [{Dia1 i2 }]
δ ψ̃n∗ (r) δ ψ̃n∗ (r)
δTs [{ψ̃n }]
=
δ ψ̃n∗ (r)
δExc [ñ] δUH [ρ̃] δñ(r′ )
Z  
+ dr′ +
δñ(r′ ) δñ(r′ ) δ ψ̃n∗ (r)
"Z #
XX a ′
′ δUH [ρ̃] δ Z̃ (r ) δ∆E a δDia1 i2
+ dr +
a i1 i2 δ Z̃ a (r′ ) δDia1 i2 δDia1 i2 δ ψ̃n∗ (r)
= fn (− 21 ∇2 )ψn (r)
Z
+ dr′ [vxc [ñ](r′ ) + uH [ρ̃](r′ )] fn δ(r − r′ )ψ̃n (r′ )
" #
XX Z
′ ′
X
a a ′ δ∆E a
+ dr uH [ρ̃](r ) ∆Li1 i2 g̃L (r ) + a fn p̃ai1 (r)Pni
a

a i i
δD i i
2
1 2 L 1 2

δExc [n]
where vxc [n](r) = is the usual local (LDA) or semilocal (GGA) exchange
δn(r)
R ′ n(r′ )
correlation potential, and uH [n](r) = δU H [n]
δn(r) = dr |r−r′ | is the usual Hartree
potential.
From these results, we can write down the transformed Hamiltonian as
XX
e = − 1 ∇2 + uH [ρ̃] + vxc [ñ] +
H
b
|p̃ai1 ⟩∆Hia1 i2 ⟨p̃ai2 |, (56)
2
a i1 i2

where the nonlocal part of the Hamiltonian is given in terms of the tensor
δ∆E a
X Z
a a a
∆Hi1 i2 = ∆Li1 i2 druH [ρ̃](r)g̃L (r) +
δDia1 i2
L
Z
X X δ∆Exc
= ∆aLi1 i2 druH [ρ̃](r)g̃L
a
(r) + ∆Tia1 i2 + ∆Cia1 i2 + 2 ∆Cia1 i2 i3 i4 Dia3 i4 + .
i i
δDia1 i2
L 3 4

(57)

Note that to justify taking the derivative with respect to D only, and not its com-
plex conjugate, the symmetry properties 48 has been used to get Dia∗ 1 i2
∆Cia1 i2 i3 i4 Dia3 i4 =
a a a
Di1 i2 ∆Ci1 i2 i3 i4 Di3 i4 .

15
1.7 Forces in PAW
In the ground state, the forces on each nuclei can be calculated directly from
dE
Fa = −
dRa ( )
∂E X ∂E d|ψ̃n ⟩
=− − + h.c.
∂Ra n ∂|ψ̃n ⟩ dRa
( ) (58)
∂E X d|ψ̃n ⟩
=− − fn ϵn ⟨ψ̃n |Ŝ + h.c.
∂Ra n
dRa
∂E X dŜ
=− a
+ fn ϵn ⟨ψ̃n | a
|ψ̃n ⟩
∂R n
dR

where h.c. denotes the hermitian conjugate. To get to the second line, the chain
rule has been applied. The third line follows from the relation
∂E e ψ̃n ⟩ = fn ϵn Ŝ|ψ̃n ⟩.
= fn H| (59)
b
∂⟨ψ̃n |

The last line of 58 is obtained from the following manipulation of the orthogo-
nality condition 53

δnm = ⟨ψ̃n |Ŝ|ψ̃m ⟩

d d⟨ψ̃n | dŜ d|ψ̃m ⟩
⇒ 0= ⟨ψ̃n |Ŝ|ψ̃m ⟩ = Ŝ|ψ̃m ⟩ + ⟨ψ̃n | |ψ̃m ⟩ + ⟨ψ̃n |Ŝ
dRa dRa dRa dRa (60)
d⟨ψ̃n | dŜ
⇔ a
Ŝ|ψ̃m ⟩ + h.c. = −⟨ψ̃n | |ψ̃m ⟩
dR dRa
From the expression for the overlap operator 54, it follows that

d|p̃ai1 ⟩ a
 
dŜ X
= ∆Sia1 i2 ⟨p̃ | + h.c. (61)
dR a
i i
dRa i2
1 2

which, when inserted in 58, gives the force expression

a
dp̃ai1 a
 
∂E a∗ dp̃i2
X X
Fa = − + f ϵ
n n ∆S a
i1 i2 P ni1 ⟨ |ψ̃ n ⟩ + ⟨ ψ̃n | ⟩P (62)
∂Ra n i i
dRa dRa ni2
1 2

In the case of standard xc approximations, the dependence of the total energy

on the nuclei coordinates is
∂E
Z ′ X ∂E ∂Dia i Z X δE ∂g̃ a (r′ )
′ δE ∂ ñ(r )
a
= dr ′ a
+ a
1 2
a
+ dr′ L
a (r′ ) ∂Ra
∂R δñ(r ) ∂R i1 i2
∂D i1 i2 ∂R δg̃L
L
a ′ a
∂Di1 i2 ∂g̃ a (r′ )
Z Z
∂ ñ (r ) X X
= dr′ [uH (r′ ) + vxc (r′ )] c a + ∆Hia1 i2 a
+ dr′ uH (r′ )QL L a
∂R i i
∂R ∂R
1 2 L
(63)

16
giving the force expression
( )
∂ ñac (r′ ) a ′
Z
a ′ ′ ′′
X ∂g̃L (r )
F = − dr [uH (r ) + vxc (r )] + uH (r ) QL
∂Ra ∂Ra
L
a
dp̃ai1 a
 
a∗ dp̃i2
X X
a a
− fn ∆Hi1 i2 − ϵn ∆Si1 i2 Pni1 ⟨ |ψ̃n ⟩ + ⟨ψ̃n | ⟩P
n i i
dRa dRa ni2
1 2

(64)

1.8 Summary
The PAW KS equation to be solved is

e ψ̃n ⟩ = ϵn Ŝ|ψ̃n ⟩
H|
b
(65)

with Ŝ, and H given by

be
XX √
Sb = 1̂ + |p̃ai1 ⟩ 4π∆a00,i1 i2 ⟨p̃ai2 | (66a)
a i1 i2
XX
H = − 21 ∇2 + uH [ρ̃](r) + vxc [ñ](r) + |p̃ai1 ⟩∆Hia1 i2 ⟨p̃ai2 | (66b)
be
a i1 i2

where
Z
X X δ∆Exc
∆Hia1 i2 = ∆aLi1 i2 a
druH [ρ̃](r)g̃L (r)+∆Tia1 i2 +∆Cia1 i2 +2 ∆Cia1 i2 i3 i4 Dia3 i4 +
i3 i4
δDia1 i2
L
(67)
The total energy can then be evaluated by
X
E = Ts [{ψ̃n }] + UH [ρ̃] + Exc [ñ] + ∆E a (68)
a

with ∆E a given by
X X
∆E a = Tca + ∆Tia1 i2 + ∆Cia1 i2 Dia1 i2 + Dia∗ ∆Cia1 i2 i3 i4 Dia3 i4 +∆Exc
a
({Dia1 i2 })


1 i2
i1 i2 i1 i2 i3 i4
(69)
Having solved the eigenvalue problem 65, the eigenvalues are known. This
can be used to determine, for example, the kinetic energy of the pseudo wave
functions, Ts [ñ], without doing the explicit (and P computationally costly) com-
putation. This can be seen by operating with n fn ⟨ψ̃n | on eq. 66b to get:
X Z XX
Ts [{ψ̃n }] = fn ϵn − dr[ñ(r)−ñc (r)] [uH [ρ̃](r) + vxc [ñ](r)]− ∆Hia1 i2 Dia1 i2
n a i1 i2
(70)

2 Implementing PAW
For an implementation of PAW, one must specify a large number of data for
each chemical element. This constitutes a data set which uniquely determines

17
how the on-site PAW transformation works, at the site of the specific atom.
For the generation of such data sets, one needs an atomic DFT program, with
which basis sets can be generated. How to perform DFT calculations efficiently
on an isolated atom will be discussed in the first section of this chapter, and
the actual choice of data set parameters will be discussed in the second. The
atomic DFT program plays the additional role of a small test program, against
which implementations in the full PAW program can be tested.

2.1 Atoms
If we consider the Kohn-Sham equation for an isolated atom, (described by a
non spin-dependent Hamiltonian), it is well known that the eigenstates can be
represented by the product

ϕiσi (rσ) = Rj (r)YL (r̂)χσi (σ) (71)

where Rj are real radial function, and YL are the (complex valued) spherical
harmonics. i = (n, l, m), j = (n, l), and L = (l, m).
Assuming identical filling of all atomic orbitals, i.e. fiσ = fj , the density
becomes
XX X 2l + 1
n(r) = fj |ϕiσi (rσ)|2 = 2 fj |Rj (r)|2 (72)
i σ j
4π
i

where the first factor of 2 comes from the sum over spin, and the second factor
from the sum over the magnetic quantum number using that
X 2l + 1
|Ylm |2 = (73)
m
4π

The identical filling of degenerate states is exact for closed shell systems,
and corresponds to a spherical averaging of the density for open shell systems.
Determining potentials in a spherical coordinate system is usually done by
exploiting the expansion of the Coulomb kernel

1 X 4π rl
= <
Y ∗ (r̂)YL (r̂′ )
l+1 L
(74)
|r − r′ | 2l + 1 r>
L

with r< = min(r, r′ ) and r> = max(r, r′ ). Using this it is seen that for any
density with a known angular dependence, e.g. the density R(r)YL (r̂), the
potential can be determined by

R(r′ )YL (r̂′ )

Z
v[R(r)YL (r̂)](r) = dr′
|r − r′ |
Z ∞ l
4π r<
= YL (r̂) r′2 dr′ R(r′ ) l+1
2l + 1 0 r>
4π
Z r  r′ l+1 Z ∞  l 
′ ′ ′ ′ ′ ′ r
= YL (r̂) dr R(r )r + dr R(r )r
2l + 1 0 r r r′
(75)

if the angular dependence is not a spherical harmonic, one can always do a

multipole expansion, and use the above expression on the individual terms.

18
 In the case of a radial density n(r) = n(r), the Hartree potential becomes
Z ∞
4π r ′
Z
uH (r) = dr n(r′ )r′2 + 4π dr′ n(r′ )r′ (76)
r 0 r

A purely radial dependent density also implies that the xc-potential is a radial
function. Using this, the entire KS equation can be reduced to a 1D problem in
r, while the angular part is treated analytically.

2.1.1 The Radial Kohn-Sham Equation

For a spherical KS potential, and using that YL are eigenstates of the Laplacian,
the KS equation can be reduced to the simpler one-dimensional second order
eigenvalue problem

1 d2
 
1 d l(l + 1)
− − + + vs (r) Rj (r) = ϵj Rj (r) (77)
2 dr2 r dr 2r2

If we introduce the radial wave function uj (r) defined by

rRj (r) = uj (r) (78)

the KS equation can be written as

 
′′ l(l + 1)
uj (r) + 2ϵj − 2vs (r) − uj (r) = 0 (79)
r2

which is easily integrated using standard techniques. See e.g. [6, chapter 6].

2.2 The Atomic Data Set of PAW

The very large degree of freedom when choosing the functions defining the PAW
transformation means that the choice varies a great deal between different im-
plementations. In any actual implementation expansions are obviously finite,
and many numerical considerations must be made when choosing these basis
sets, to ensure fast and reliable convergence. This section provides an overview
of the information needed for uniquely defining the PAW transformation, and
the level of freedom when choosing the parameters.

The Partial Waves

The basis functions, ϕai , for the expansion of |ψn ⟩ should be chosen to ensure
a fast convergence to the KS wave function. For this reason we choose the
partial waves as the eigenstates of the Schrödinger equation for the isolated
spin-saturated atoms. Thus the index i is a combination of main-, angular-,
and magnetic quantum number, (n, l, m). And the explicit form is

ϕai (r) = ϕanl (r)Ylm (r̂)

where ϕanl (r) are the solutions of the radial KS Schrödinger equation 77, and
Ylm are the spherical harmonics. For convenience we choose ϕai (r) to be real, i.e.
we use the real spherical harmonics instead of the complex valued. This choice
of partial waves implies that the smooth partial waves and the smooth projector

19
functions can also be chosen real, and as products of some radial functions and
the same real spherical harmonic.
Note that including unbound states of the radial KS equation in the partial
waves is not a problem, as the diverging tail is exactly canceled by the smooth
partial waves. In practice we only integrate the KS equation outward from the
origin to the cutoff radius for unbound states, thus making the energies free
parameters. In principle the same could be done for the bound states, but in
gpaw, the energies of bound states are fixed by making the inward integration
for these states and doing the usual matching (see e.g. [6, chapter 6]), i.e. the
energies are chosen as the eigen energies of the system.

The Smooth Partial Waves

The smooth partial waves ψ̃ia (r) are per construction identical to the partial
waves outside the augmentation sphere. Inside the spheres, we can choose them
as any smooth continuation. Presently gpaw uses simple 6’th order polynomials
of even powers only (as odd powers in r results in a kink in the functions at
the origin, i.e. that the first derivatives are not defined at this point), where
the coefficients are used to match the partial waves smoothly at r = rc . Other
codes uses Bessel functions[4] or Gaussians.

The Smooth Projector Functions

The smooth projector functions are a bit more tricky. Making them orthonormal
to ϕ̃ai (r) is a simple task of applying an orthonormalization procedure. This is
the only formal requirement, but in any actual implementation all expansions
are necessarily finite, and we therefore want them to converge as fast as possible,
so only a few terms needs to be evaluated.
A popular choice is to determine the smooth projector functions according
to
|p̃ai ⟩ = − 12 ∇2 + ṽs − ϵi |ϕ̃ai ⟩


(80)
or equivalently

1 d2
 
1 d l(l + 1)
p̃aj (r) = − − + + ṽs (r) − ϵj ϕ̃aj (r) (81)
2 dr2 r dr 2r2

where ṽs (r) is the smooth KS potential uH [ρ̃](r) + vxc [ñ](r). And then enforce
the complementary orthogonality condition ⟨p̃aj |ϕ̃aj′ ⟩ = δj,j ′ inside the augmenta-
tion sphere, e.g. by a standard Gram-Schmidt procedure. Using this procedure
ensures that the reference atom is described correctly despite the finite number
of projectors.

The Smooth Compensation Charge Expansion Functions

a
The smooth compensation charges g̃L (r), are products of spherical harmonics,
a
and radial functions g̃l (r) satisfying that
Z
drrl YL (r̂)g̃L
a
′ (r) = δLL′ (82)

20
In gpaw the radial functions are chosen as generalized Gaussian according to
(here shown for Ra = 0):

1 l! a 2
a
g̃L (r) = g̃la (r)YL (r̂) , g̃la (r) = √ (4αa )l+3/2 rl e−α r (83)
4π (2l + 1)!
where the atom-dependent decay factor α is chosen such that the charges are
localized within the augmentation sphere.

The Core- and Smooth Core Densities

The core density follows directly from the all electron partial waves by
core
X core
X
nc (r) = |ϕi (r)|2 = 2(2l + 1)|ϕj (r)|2 /4π (84)
i j

The smooth core densities ñac (r) are like the smooth partial waves expanded
in a few (two or three) Bessel functions, Gaussians, polynomials or otherwise,
fitted such that it matches the true core density smoothly at the cut-off radius.

The Localized Potential

An additional freedom in PAW is that for any operator L,
b localized within the
augmentation spheres, we can exploit the identity 8
X
|ϕ̃ai ⟩⟨p̃ai | = 1 (85)
i

valid within the spheres, to get

XX
L
b= |p̃ai1 ⟩⟨ϕ̃ai1 |L|
b ϕ̃ai ⟩⟨p̃ai |
2 2
a i1 i2

so for any potential v̄(r) = a v̄ a (r − Ra ) localized within the augmentation

P
spheres (i.e. v̄ a (r) = 0 for r > rca ) we get the identity
Z X XZ
a
0 = drñ(r) v̄ (r) − drña v̄ a (r)
a a

This expression can be used as an ‘intelligent zero’. Using this, we can make
the replacement of the smooth potential

ṽeff (r) = uH [ρ̃](r) + vxc [ñ](r) → ṽeff (r) = uH [ρ̃](r) + vxc [ñ](r) + v̄(r) (86)

if we at the same time add

X
Ba + Bia1 i2 Dia1 i2 (87)
i1 i2

to the energy corrections ∆E a , where

Z Z
B a = − drñac v̄ a (r) and ∆Bia1 i2 = − drϕ̃ai1 ϕ̃ai2 v̄ a (r) (88)

21
This also implies that Bia1 i2 should be added to ∆Hia1 i2 .
The advantage of doing this is that the Hartree potential and the xc-potential
might not be optimally smooth close to the nuclei, but if we define the localized
potential properly, the sum of the three potentials might still be smooth. Thus
one can initially evaluate uH [ρ̃](r) and vxc [ñ](r) on an extra fine grid, add v̄(r)
and then restrict the total potential to the coarse grid again before solving the
KS equation.
The typical way of constructing the localized potentials v̄ a is by expanding
it in some basis, and then choosing the coefficients such that the potential
uH [ρ̃](r) + vxc [ñ](r) + v̄(r) is optimally smooth at the core for the reference
system.
Inclusion of v̄ a (r) changes the forces on each atom only through the redefi-
nitions of ṽeff (r) and ∆Hia1 i2 .

Summary
When constructing a data set for a specific atom, one must specify the following
quantities, all defined within the augmentation spheres only:
1. ϕai from radial KS equation
2. ϕ̃ai by appropriate smooth continuation of ϕai
3. p̃ai from equation 80
′
a
drrl g̃L
a
R
4. g̃L localized within r < rc , and satisfying (r)YL′ (r −
d Ra ) = δLL′
5. nac follows from ϕai by 84

6. ñac by appropriate smooth continuation of nac

7. v̄ a localized within r < rca , otherwise freely chosen to make ṽeff optimally
smooth for the reference system
The adjustable parameters besides the shape of ϕ̃a , g̃L
a
, v̄ a , and ñac are

1. Cut-off radii rca (which can also depend on i)

2. Frozen core states
3. Number of basis set functions (range of index i on ϕai , ϕ̃ai , and p̃ai )
4. Energies of unbound partial waves

Choosing these parameters in such a way that the basis is sufficient for the
description of all possible environments for the specific chemical element, while
still ensuring a smooth pseudo wave function is a delicate procedure, although
the optimal parameter choice is more stable than for e.g. norm conserving or
ultra soft pseudopotentials.
Once the quantities above have been constructed, all other ingredients of
the PAW transformation follows, such as ∆a , ∆aLii′ , Tca , ∆Tia1 i2 , ∆C a , ∆Cia1 i2 ,
∆Cia1 i2 i3 i4 , B a , and Bia1 i2 . The first two are needed for the construction of the
compensation charges and the overlap operator, and the rest for determining
the Hamiltonian, and for evaluating the total energy.

22
3 Non-standard Quantities
The preceding sections have described the details of making a standard DFT
scheme work within the PAW formalism. This section will focus on what the
PAW transform does to quantities needed for post-processing or expansions to
DFT.
It is a big advantage of the PAW method, that it is formally exactly equiv-
alent to all-electron methods (with a frozen core) but is computationally com-
parable to doing pseudopotential calculations. In pseudopotential approaches,
projecting out the core region is handled by a static projection kernel, while in
PAW this projection kernel is dynamically updated during the SCF-cycle via an
expansion of the core region in a local atomic basis set. This has the drawback
for the end user, the equations for all quantities most be modified to account
the dual basis set description.

3.1 The External Potential in PAW

As an example of the principle in accommodating expressions to the PAW for-
malism, we will here consider the application of an external potential in DFT.
Without the PAW transformation, this addition is trivial, as the desired
potential, vext (r), should simply be added to the effective KS potential, and
the total
R energy adjusted by the energy associated with the external potential
Eext = drn(r)vext (r).
In PAW, the density decomposes into pseudo and atomic parts, so that
Z XZ
Eext = drñ(r)vext (r) + dr [na (r) − ña (r)] vext (r).
a
R
Implying that both a pseudoR energy contribution Ẽext = drñ(r)vext (r) and
a
atomic corrections ∆Eext = dr [na (r) − ña (r)] vext (r) should be added to the
total energy.
In PAW, the Hamiltonian has the structure:
1 ∂E XX
H= = H̃ + |p̃ai1 ⟩∆Hia1 i2 ⟨p̃ai2 |
fn |ψn ⟩ ∂⟨ψn | a i i 1 2

In our case, the extra contributions due to the external potential are:

H̃ext (r) = vext (r)

and Z n o
∆Hia,ext
1 i2
= drvext (r) ϕai1 (r)ϕai2 (r) − ϕ̃ai1 (r)ϕ̃ai2 (r) (89)

Thus we can write the atomic energy contribution as:

Z " #
X n o
a a a a a a a a
∆Eext = drvext (r) nc (r) − ñc (r) + Di1 i2 ϕi1 (r)ϕi2 (r) − ϕ̃i1 (r)ϕ̃i2 (r)
i1 i2
Z X
= drvext (r) [nac (r) − ñac (r)] + Dia1 i2 ∆Hia,ext
1 i2
i1 i2

23
To evaluate the first term in the last line, the external potential should be
a
expanded in some radial function at each nuclei e.g. the gaussians g̃L , as the
integral of these with the core densities is already precalculated.
For example, a zero-order (monopole) expansion, equivalent to the assump-
tion
vext (r) ≈ vext (Ra ) , for |r − Ra | < rca
Leads to the expression:
a
√
∆Eext = vext (Ra )( 4πQa00 + Z a )
√
∆Hia,ext
1 i2
= vext (Ra ) 4π∆a00,i1 i2

Linear external potentials (corresponding to a homogeneous applied electric

field) can be handled exactly by doing an expansion to first order. This has
been used in gpaw in e.g. the paper [7].

3.2 All-electron Density

During the self-consistency cycle of DFT, the all-electron quantities are at all
times available in principle. In practise, they are never handled directly, but
rather in the composite basis representation of a global pseudo description aug-
mented by local atomic basis functions. For some post processing purposes it
can however be desirable to reconstruct all-electron quantities on a single regular
grid.
As an example, consider the all-electron density, which can formally be re-
constructed by
" #
X X  
a a a a a a
n(r) = ñ(r) + nc (r) + Di1 i2 ϕi1 (r)ϕi2 (r) − ϕ̃i1 (r)ϕ̃i2 (r) .
a i1 i2

Transferring this to a uniform grid will of coarse re-introduce the problem

of describing sharply peaked atomic orbitals on a uniform grid, but as it is only
needed for post processing, and not in the self-consistency, it can be afforded to
interpolating the pseudo density to an extra fine grid, before adding the summed
atomic corrections from the radial grid.
One common use of the all-electron density is for the application of Bader
analysis[8]. The advantage of applying this to the all-electron density instead
of the pseudo density, is that it can be proved that the total electron density
only has maxima’s at the nuclei, such that there will only be one Bader volume
associated with each atom. This does not hold for the pseudo density, which can
result in detached Bader volumes. In addition, the dividing surfaces found if
applied to the pseudo density will be wrong if these intersect the augmentation
sphere.
In practice, the reconstructed total density will not integrate correctly due
to the inaccurate description of a uniform grid in the core regions of especially
heavy elements. But as the numerically exact value of the integral P over the
a
atomic corrections are known from the radial grid description (= 4π ij Dij ∆aL,ij
), this deficiency can easily be remedied. As long as the density is correctly de-
scribed at the dividing surfaces, these will still be determined correctly.

24
3.3 Wannier Orbitals
When constructing Wannier functions, the only quantities that needs to be sup-
plied by the DFT calculator are the integrals ZnG1 n2 = ⟨ψn1 |e−G·r |ψn2 ⟩, where
G is one of at most 6 possible (3 in an orthorhombic cell) vectors connecting
nearest neighbor cells in the reciprocal lattice.[9, 10]
When introducing the PAW transformation, this quantity can be expressed
as
XX  
ZnG1 n2 = ⟨ψ̃n1 |e−G·r |ψ̃n2 ⟩+ Pna∗1 i1 Pna2 i2 ⟨ϕai1 |e−G·r |ϕai2 ⟩ − ⟨ϕ̃ai1 |e−G·r |ϕ̃ai2 ⟩ .
a i1 i2

Even for small systems, the phase of the exponential of the last integral does not
vary significantly over the region of space, where p̃ai is non-zero. The integral in
the last term can therefore safely be approximated by
a X √
e−G·R Pna∗1 i1 Pna2 i2 4π∆a00,i1 i2 .
i1 i2

3.4 Local Properties

This section describes quantities that can somehow be related to a specific atom.
As the PAW transform utilizes an inherent partitioning of space into atomic
regions, such quantities are usually readily extractable from already determined
a
atomic attributes, such as the density matrices or the projector overlaps Pni ,
which are by construction simultaneous expansion coefficients of both the pseudo
and the all-electron wave functions inside the augmentation spheres.

3.4.1 Projected Density of States

The projection of the all electron wave functions onto the all electron partial
waves, (i.e. the all electron wave functions of the isolated atoms) ϕai , is within
the PAW formalism given by
XX ′
 ′ ′ ′ ′
 ′
⟨ϕai |ψn ⟩ = ⟨ϕ̃ai |ψ̃n ⟩ + ⟨ϕ̃ai |p̃ai1 ⟩ ⟨ϕai1 |ϕai2 ⟩ − ⟨ϕ̃ai1 |ϕ̃ai2 ⟩ ⟨p̃ai2 |ψ̃n ⟩ (90)
a ′ i1 i2

Using that projectors and pseudo partial waves form a complete basis within
the augmentation spheres, this can be re-expressed as
′ ′
a′
XX
⟨ϕai |ψn ⟩ = Pni
a
+ ⟨ϕ̃ai |p̃ai1 ⟩∆Sia1 i2 Pni 2
(91)
a′ ̸=a i1 i2
′
if the chosen orbital index ‘i‘ correspond to a bound state, the overlaps ⟨ϕ̃ai |p̃ai1 ⟩,
a′ ̸= a will be small, and we see that we can approximate
⟨ϕai |ψn ⟩ ≈ ⟨p̃ai |ψ̃n ⟩ (92)
a
The coefficients Pni
= ⟨p̃ai |ψ̃n ⟩,
can thus be used as a qualitative measure of the
local character of the true all electron wave functions. As the coefficients are
already calculated and used in the SCF cycle, it involves no extra computational
cost to determine quantities related directly to these.
These can be used to define an atomic orbital projected density of states
X
a 2
ni (ε) = δ(ε − ϵn ) |Pni | . (93)
n

25
3.4.2 Local Magnetic Moments
As the projection coefficients are simultaneous expansion coefficients of the
pseudo and the all-electron wave functions inside the augmentation spheres, it
can be seen that inside these, the all-electron density is given by (for a complete
set of partial waves)
X
n(r) = Dia1 i2 ϕai1 (r)ϕai2 (r) + nac (r) , |r − Ra | < rca . (94)
i1 i2

This can be used to assign a local magnetic moment to each atom according
to X
Ma = ∆Nia1 i2 Dia1 i2 (↑) − Dia1 i2 (↓) ,
 
i1 i2

where ∆N is an integration over products of AE waves truncated to the interior

of the augmentation sphere
Z
∆Nia1 i2 = drϕai1 (r)ϕai2 (r).
r∈Ωa
R
Note that this will not add up to the total magnetic moment dr(n↑ (r) −
n↓ (r)), due to the interstitial space between augmentation spheres, and must
be scaled if this is desired.

3.4.3 LDA + U
The atom projected density matrix Dia1 i2 can also be used to do LDA + U
calculations. The gpaw implementation follows the LDA + U implementation
in VASP[11], which is based on the particular branch of LDA + U suggested by
Dudarev et al.[12], where you set the effective (U-J) parameter. The key notion
is that from 94 one can define an (valence-) orbital density matrix

ρ̂ai1 i2 = |ϕai1 ⟩Dia1 i2 ⟨ϕai2 |.

Thus doing LDA + U is a simple matter of picking out the d-type elements of
Da , and adding to the total energy the contribution
d type
!
UX X X
Dia1 i1 − Dia1 i2 Dia2 i1 (95)
2 a i i
1 2

and adding the gradient of this, wrt. Dia1 i2 , to the atomic Hamiltonian

UX a
|p̃i1 ⟩ δi1 i2 − 2Dia1 i2 ⟨p̃ai2 |


(96)
2 a

3.5 Coulomb Integrals

When trying to describe electron interactions beyond the level of standard (semi-
) local density approximations, one will often need Coulomb matrix elements of
the type
drdr′ ∗
ZZ
Knn′ ,mm′ = (nnn′ |nmm′ ) := n ′ (r)nmm′ (r′ ), (97)
|r − r′ | nn

26
where the orbital pair density nnn′ (r) = ψn∗ (r)ψn′ (r).
Such elements are needed in some formulations of vdW functionals (although
not the one implemented in gpaw), in linear-response TDDFT (see e.g. [13])
where only pair densities corresponding to electron-hole pairs are needed, in
exact exchange or hybrid functionals (see next section) where only elements of
the form Knn′ ,nn′ where both indices correspond to occupied states, are needed,
and for GW calculations (see e.g. [14]), where all elements are needed.
Introducing the PAW transformation in 97, the pair densities partition ac-
cording to X
nnn′ (r) = ñnn′ (r) + (nann′ (r) − ñann′ (r)) (98)
a

with the obvious definitions

X X
ñnn′ =ψ̃n∗ ψ̃n′ nann′ = a∗ a
Pni P ′ ϕa∗ ϕa∗
1 n i2 i1 i2 ñann′ = a∗ a
Pni P ′ ϕ̃a∗ ϕ̃a∗
1 n i2 i1 i2 . (99)
i1 i2 i1 i2

Exactly like with the Hartree potential, direct insertion of this in 97 would, due
to the non-local nature of the Coulomb kernel, lead to undesired cross terms
between different augmentation spheres. As before, such terms can be avoided
a
by introducing some compensation charges, Z̃nn ′ , chosen such that the potential
a a a
of nnn′ − ñnn′ − Z̃nn′ are zero outside their respective augmentation spheres.
This is achieved by doing a multipole expansion and requiring the expansion
coefficients to be zero, and entails a compensation of the form
X X
a
Z̃nn′ (r) = QaL,nn′ g̃L
a
(r), QaL,nn′ = ∆aL,i1 i2 Pni
a∗ a
P ′
1 n i2
(100)
L i1 i2

(the constants ∆aL,i1 i2 are identical to those in 41b).

Introduction of such compensation charges makes it possible to obtain the
clean partitioning
X X
a
Knn′ ,mm′ = (ρ̃nn′ |ρ̃mm′ ) + 2 Pmi P a∗ ∆Cia1 i2 i3 i4 Pna∗′ i3 Pm
1 ni2
a
′i .
4
(101)
a i1 i2 i3 i4

a
Here the last term is a trivial functional of the expansion coefficients Pni in-
a
volving only the constants ∆Ci1 i2 i3 i4 already precalculated for the atomic cor-
rections to the Coulomb energy 47. The only computationally demanding term
relates to thePCoulomb matrix element of the smooth compensated pair densities
a
ρ̃ij = ñij + a Z̃ij , which are expressible on coarse grids.
The formally exact partitioning 101 makes it possible, at moderate computa-
tional effort, to obtain Coulomb matrix elements in a representation approaching
the infinite basis set limit. In standard implementations, such elements are usu-
ally only available in atomic basis sets, where the convergence of the basis is
problematic. At the same time, all information on the nodal structure of the
all-electron wave functions in the core region is retained, which is important
due the non-local probing of the Coulomb operator. In standard pseudopoten-
tial schemes, this information is lost, leading to an uncontrolled approximation
to Knn′ ,mm′ .
As a technical issue, we note that integration over the the Coulomb kernel
1/|r − r′ | is done by solving the associated Poisson equation, as for the Hartree
potential, whereby the calculation of each element can be efficiently parallelized

27
 R
using domain decomposition. The integral drρ̃nn′ (r) = δnn′ shows that the
compensated pair densities ρ̃nn have a non-zero total charge, which is problem-
atic for the determination of the associated potential. For periodic systems,
charge neutrality is enforced by subtracting a homogeneous background charge,
and the error so introduced is removed to leading order (V −1/3 where V the
the volume of the simulation box) by adding the potential of a missing probe
charge in an otherwise periodically repeated array of probe charges embedded
in a compensating homogeneous background charge. This can be determined
using the standard Ewald technique, and corresponds to a rigid shift of the
potential. For non-periodic systems, all charge is localized in the box, and the
Poisson equation can be solved by adjusting the boundary values according to a
multipole expansion of the pair density with respect to the center of the simula-
tion box. A monopole correction is correct to the same order as the correction
for periodic cells, but the prefactor on the error is much smaller, and leads to
converged potentials even for small cells.

3.5.1 Exact Exchange

The EXX energy functional is given by
1X
Exx == − fn fm δσn ,σm Knm,nm . (102)
2 nm

Terms where n and m both refer to valence states transform in PAW as in

equation 101. Terms where either index refers to a core orbital can be reduced
a
to trivial functionals of Pni , resulting in (see e.g. [15])
val
1X
Exx = − fn fm δσn ,σm (ρ̃nm |ρ̃nm )
2 nm
" #
X X X X
− Dia∗
1 i3
(σ)∆Cia1 i2 i3 i4 Dia2 i4 (σ) + Dia1 i2 Xia1 i2 + Exx
a,c-c
.
a σ i1 i2 i3 i4 i1 i2
(103)
The term involving the ∆C a tensor is the PAW correction for the valence-valence
interaction, and is similar to the correction in the equivalent expression for the
Hartree energy, except that the order of the indices on the density matrices
are interchanged. The term involving the X a tensor represents the valence-core
a,c-c
exchange interaction energy. Exx is simply the (constant) exchange energy of
the core electrons.
The system independent Hermitian tensor Xia1 i2 is given by:
core Z Z
1X ϕai1 (r)ϕa,core
α (r)ϕai2 (r′ )ϕa,core
α (r′ )
Xia1 i2 = drdr′
2 α |r − r′ |
core X
! ZZ
l
X 4π X r<
= L L
GL1 Lc GL2 Lc drdr′ l+1 uaj1 (r)uajc (r)uajc (r′ )uaj2 (r′ ).
jc
2l + 1 mmc r >
l
(104)
Although the valence-core interaction is computationally trivial to include,
it is not unimportant, giving rise to shifts in the valence eigenvalues of up

28
to 1eV (though only a few kcal/mol in atomization energies), and we note
that this contribution is unavailable in pseudopotential schemes. The core-
core exchange is simply a reference energy, and will not affect self-consistency
or energy differences.
For the iterative minimization schemes used in real-space and plane wave
codes, the explicit form of the non-local Fock operator v NL (r, r′ ) is never needed,
and would indeed be impossible to represent on any realistic grid. Instead
only the action of the operator on a state is needed. As with the Hamiltonian
operator, the action on the pseudo waves is derived via the relation fn v̂ NL |ψ̃n ⟩ =
∂Exx /∂⟨ψ̃n |. Referring to [15] for a derivation, we merely state the result
X
v̂ NL |ψ̃n ⟩ = fm ṽnm (r)|ψ̃m ⟩
m
" ! #
XX X X
+ |p̃ai1 ⟩ a
vnm,i P a − Xia1 i2 Pni
1 i2 mi2
a
2
−2 Cia1 i3 i2 i4 Dia3 i4 a
Pni 2
a i1 i2 m i3 i4
(105)

where ṽnm is the solution of ∇2 ṽnm (r) = −4π ρ̃nm (r), and vnm,i
a a a
P R
1 i2
= L ∆Li1 i2 drg̃L (r)ṽnm (r).
Again the computationally demanding first term is related to smooth pseudo
quantities only, which can be accurately represented on coarse grids, making it
possible to do basis set converged self-consistent EXX calculations at a relatively
modest cost. Applying the Fock operator is however still expensive, as a Poisson
equation must be solved for all pairs of orbitals.
As a technical consideration, note that the effect of the atomic corrections
due to valence-valence, valence-core, and core-core exchange interactions can
simply be incorporated into the standard equations by redefining equations 47,
46, and 45 respectively, which will also take care of the last two terms in the
Fock operator above. The introduction of the pair orbital compensation charges
does however lead to a non-trivial correction to the Fock operato; the term
a
proportional to vnm,i 1 i2
. This term also leads to a distinct contribution when
calculating the kinetic energy via the eigenvalues as done in equation 70. The
additional term (besides those related to redefining 45–47)
" Z #
X XX
fn fm δσn ,σm drṽnm (r)ψ̃n∗ (r)ψ̃m (r) − a
Pni P a va
1 mi2 nm,i1 i2
, (106)
nm a i1 i2

should be added to the right hand side of 70 on inclusion of exact exchange.

In a similar fashion, the compensation charges leads to an additional force
contribution in equation 64 given by
(Z
a ′
X X X ∂g̃L (r )
a
Fxx = fn fnm δσn σm dr′ ṽnm (r′ ) a∗ a
Pni P
1 mi2
∆ Li i
1 2
nm i1 i2
∂Ra
L
a )
dp̃ai1 a

a∗ dp̃i2
X
a
+ vn1 n2 i1 i2 Pni1 ⟨ |ψ̃m ⟩ + ⟨ψ̃n | ⟩P .
i i
dRa dRa mi2
1 2

(107)

29
3.5.2 Optimized Effective Potential
The optimized effective potential (OEP) method, is a way of converting the
non-local Fock operator v̂xNL into a local form v̂xL = vxL (r).
One way to derive the OEP equations in standard KS-DFT, is to use pertur-
bation theory along the adiabatic connection (Görling-Levy perturbation theory
[16]).
On converting the OEP equation to the PAW formalism, it should be re-
membered that local potentials in PAW transform to a local pseudo part plus
non-local atomic corrections. Hence we want to arrive at a potential of the form
XX
v̂xL = ṽxL (r) + |p̃ai1 ⟩∆via1 i2 ⟨p̃ai2 |, (108)
a i1 i2

where both the pseudo part ṽxL as well as the coefficients ∆via1 i2 should be
determined.
The derivation is more or less straight forward, if one remembers the the
PAW KS equation is a generalized eigenvalue problem, that the variational
quantity is the pseudo orbitals, and that the first order shift in the density has
both a pseudo and an atomic part. The result is
X X ⟨ψ̃m |v̂xNL − v̂xL |ψ̃n ⟩
fn ψ̃n∗ (r) ψ̃m (r) + c.c. = 0 (109a)
n
ϵn − ϵm
m̸=n
X
a∗
X
a ⟨ψ̃m |v̂xNL − v̂xL |ψ̃n ⟩
fn Pni Pmi + c.c. = 0 (109b)
n
1 2
ϵn − ϵm
m̸=n

where v̂xNL is the non-local exchange operator of equation 105 and v̂xL is the local
version in 108.
These can be solved iteratively starting from a local density-function ap-
proximation to the exchange potential in the spirit of [17].
It might seem that OEP is just extra work on top of the already expensive
non-local operator, but it can in some cases be faster, as the number of SCF
iterations in the KS cycle are greatly reduced.

30
References
[1] P. E. Blöchl. Projector augmented-wave method. Physical Review B,
50(24):17953–17979, Dec 1994.
[2] The real-space PAW-DFT code GPAW is part of the CAMP
Open-Source (CAMPOS) project. GPAW is freely available at
https://wiki.fysik.dtu.dk/gpaw.
[3] P. E. Blöchl, C. J. Först, and J. Schimpl. Projector augmented wave
method: ab-initio molecular dynamics with full wave functions. Bulletin of
Materials Science, 26:33–41, 2003.

[4] G. Kresse and D. Joubert. From ultrasoft pseudopotentials to the projector

augmented-wave method. Physical Review B, 59:1758–1775, Jan 1999.
[5] M. Marsman and G. Kresse. Relaxed core projector-augmented-wave
method. The Journal of Chemical Physics, 125(10):104101, 2006.

[6] C. Fiolhais, F. Nogueira, and M. Margues, editors. A Primer in Density

Functional Theory, volume 620 of Lecture Notes in Physics. Springer, 2003.
[7] F. Yin, J. Akola, P. Koskinen, M. Manninen, and R. E. Palmer. Bright
beaches of nanoscale potassium islands on graphite in stm imaging. Physical
Review Letters, 102(10):106102, 2009.

[8] W Tang, E Sanville, and G Henkelman. A grid-based bader analysis

algorithm without lattice bias. Journal of Physics: Condensed Matter,
21(8):084204 (7pp), 2009.
[9] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen. Partly occupied
wannier functions: Construction and applications. Physical Review B,
72(12):125119, September 2005.
[10] A Ferretti, A Calzolari, B Bonferroni, and R Di Felice. Maximally localized
wannier functions constructed from projector-augmented waves or ultrasoft
pseudopotentials. Journal of Physics: Condensed Matter, 19(3):036215
(16pp), 2007.
[11] A. Rohrbach, J. Hafner, and G. Kresse. Molecular adsorption on the surface
of strongly correlated transition-metal oxides: A case study for co/nio(100).
Physical Review B, 69(7):075413, Feb 2004.
[12] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P.
Sutton. Electron-energy-loss spectra and the structural stability of nickel
oxide: An lsda+u study. Phys. Rev. B, 57(3):1505–1509, Jan 1998.
[13] Michael Walter, Hannu Häkkinen, Lauri Lehtovaara, Martti Puska, Jussi
Enkovaara, Carsten Rostgaard, and Jens Jørgen Mortensen. Time-
dependent density-functional theory in the projector augmented-wave
method. The Journal of Chemical Physics, 128(24):244101, 2008.
[14] C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen. Fully self-consistent
gw calculations for molecules. Phys. Rev. B, 81(8):085103, Feb 2010.

31
[15] Joachim Paier, Robin Hirschl, Martijn Marsman, and Georg Kresse. The
perdew–burke–ernzerhof exchange-correlation functional applied to the g2-
1 test set using a plane-wave basis set. The Journal of Chemical Physics,
122(23):234102, 2005.
[16] A. Görling and M. Levy. Exact kohn-sham scheme based on perturbation
theory. 50:196–204, 1994.
[17] S. Kümmel and J. P. Perdew. Simple iterative construction of the opti-
mized effective potential for orbital functionals, including exact exchange.
90:043004, 2003.

32