Alkauskas A., Deak P., Neugebauer J., Pasquarello A., Van de Walle Ch.G. (Eds.) Advanced Calculations for Defects in Materials: Electronic Structure Methods
Подождите немного. Документ загружается.
with caution since
^
T
nuc
and
^
1E
el
n
are non-commuting operators (both depend on the
ionic coordinates). Therefore, we have to apply the so-called Zassenhaus formula [24]:
e
bð
^
T
nuc
þ
^
1E
el
n
Þ
¼ e
b
^
T
nuc
e
b
^
1E
el
n
e
b
2
/2½
^
T
nuc
;
^
1E
el
n
e
b
3
/6ð2½
^
1E
el
n
;½
^
T
nuc
;
^
1E
el
n
þ½
^
T
nuc
;½
^
T
nuc
;
^
1E
el
n
Þ
…:
ð15:20Þ
Here, the dots denote exponentials corresponding to higher orders in b and with
increasingly nested commutators. An explicit formula for the higher order terms is
given in Ref. [24]. In Ref. [23], it is shown that exponentials corresponding to orders b
2
and higher will be small if m
e
M (m
e
¼electron mass, M ¼nucleus mass), i.e.,
under the same condition as assumed in the standard Born–Oppenheimer approx-
imation. It is therefore justified to approximate
e
bð
^
T
nuc
þ
^
1E
el
n
Þ
e
b
^
T
nuc
e
b
^
1E
el
n
; ð15:21Þ
in any case in which the Born–Oppenheimer approximation is justified. Using
Eq. (15.21), the fact that ðe
b
^
T
nuc
e
b
^
1E
el
n
Þ corresponds to a block diagonal matrix and
the invarianceproperty of the trace (which allows to choose the same basis for different
m, e.g., n ¼ v
0
; with a fixed v
0
;) yields:
Z ¼
X
v;m
L
v;m
e
b
^
T
nuc
e
b
^
1E
e1
v
L
v;m
DE
¼
X
v;m
L
v
0
;m
e
b
^
T
nuc
e
b
^
1E
e1
v
L
v
0
;m
DE
¼
X
m
L
v
0
;m
e
b
^
T
nuc
X
m
e
b
^
1E
e1
v
!
L
v
0
;m
*+
¼
X
m
L
v
0
;m
e
b
^
T
nuc
e
b
^
1F
e1
L
v
0
;m
DE
;
ð15:22Þ
with the electronic free energy defined by:
F
el
:¼k
B
T ln
X
n
exp½E
el
n
/ðk
B
TÞ: ð 15:23Þ
In order to recombine the exponentials again, we need to apply the Baker–Campbell–
Hausdorff formula which reads [24]:
e
b
^
T
nuc
e
b
^
1F
el
¼ e
bð
^
T
nuc
þ
^
1F
el
Þ
b
2
2
½
^
T
nuc
;
^
1F
el
b
3
6
½½
^
T
nuc
;
^
1F
el
;
^
1F
el
^
T
nuc
…
: ð15:24Þ
The terms in the exponential of order b
2
and higher correspond again to terms
which are small if m
e
M, thus justifying the following approximation:
e
b
^
T
nuc
e
b
^
1F
el
e
bð
^
T
nuc
þ
^
1F
el
Þ
: ð15:25Þ
In fact, as shown in Ref. [23] the terms neglected in Eq. (15.25) lead to contributions
which are of the same order as the contributions from the neglected terms in
268
j
15 Formation Energies of Point Defects at Finite Temperatures
Eq. (15.21) having however the opposite sign. Therefore, the approximations
performed in Eqs. (15.21) and (15.25) partially compensate each other.
Inserting Eq. (15.25) into Eq. (15.22) yields for the partition function and thus the
free energy:
F ¼k
B
T ln Z; ð15:26Þ
with
Z ¼
X
m
hL
m
je
bð
^
T
nuc
þ
^
1F
el
Þ
jL
m
i: ð15:27Þ
We can rewrite this to a more convenient notation as
F ¼k
B
T ln
X
m
e
b
~
E
nuc
m
; ð15:28Þ
with
^
T
nuc
þ
^
1F
el
ðfR
I
g; V; TÞ
~
L
m
¼
~
E
m
nuc
~
L
m
; ð15:29Þ
wherewehavedefined an effective nuclei Schr
€
odinger equation with eigenfunctions
~
L
m
and eigenvalues
~
E
m
nuc
. In Eq. (15.29), we have also explicitly written the dependence of F
el
on the set of nuclei coordinates fR
I
gand the crystal volume V which is a consequenc e of
the fact that each of the electronic potential energy surfaces E
el
n
depends on fR
I
gand V.
The key equations summarizing the preceding derivation are collected in Figure 15.2.
The central step which allows for a separation into an electronic, quasiharmonic,
and anharmonic part is a Taylor expansion of F
el
in Eq. (15.28) in the fR
I
garound the
T ¼ 0 K equilibrium positions fR
0
I
g:
F
el
¼ F
el
0
þ
1
2
X
k;l
u
k
u
l
q
2
F
el
qR
k
qR
l
fR
0
I
g
þOðu
3
Þ: ð15:30Þ
Here, the zeroth order term is abbreviated as F
el
0
ðV; TÞ :¼ F
el
ðfR
0
I
g; V; TÞ, k and l
run over all nuclei of the system and additionally over the three spatial dimensions for
each nucleus, and u
k
¼ R
k
R
0
k
is the displacement out of equilibrium. The expan-
sion does not contain a first order term, since such a term relates to the atomic forces
that are absent in an equilibrium structure. Each of the other terms in Eq. (15.30)
corresponds to a different excitation mechanism as emphasized in Figure 15.3. The
exact derivations are given in the following sections.
15.2.2.2 Electronic Excitations
It is convenient to decompose the electronic part F
el
0
ðV; TÞ in Eq. (15.30) into
a temperature independent part E
el
g;0
ðVÞ, which corresponds to the ground state of
the potential energy surfaces E
el
n
at fR
0
I
g, and a remainder
~
F
0
el
ðV; TÞ carrying the
temperature dependence of the electronic system, i.e., the electronic excitations:
F
el
0
ðV; TÞ¼E
el
g;0
ðVÞþ
~
F
0
el
ðV; TÞ. The reason for this separation is that
~
F
0
el
ðV; TÞ can
be accurately described with low order polynomials, while E
el
g;0
ðVÞ can be parame-
15.2 Methodology
j
269
Figure 15.3 (online colour at: www.pss-b.com) Taylor expansion of the electronic free energy F
el
as in Eq. (15.20) (see text for definitions). The color coding
emphasizes the connection with the equations in Figure 15.2. The under braces indicate the type of contribution arising from each term and the section covering the
issue.
270
j
15 Formation Energies of Point Defects at Finite Temperatures
trized using standard equations-of-state. To calculate the latter a standard DFT
approach is sufficient. For the calculation of
~
F
0
el
, one needs instead to employ the
finite temperature extension of DFT as originally developed by Mermin [25]. This
approach is implemented in typical DFTcodes and it amounts to using a Fermi–Dirac
occupation distribution for the Kohn–Sham electronic energy levels,
~
F
0
el
ðV; TÞk
B
T/2
X
i
½f
i
lnf
i
þð1f
i
Þln ð1f
i
Þ; ð15:31Þ
with self consistently determined Kohn–Sham occupation numbers f
i
¼ f
i
ðV; TÞ.
Note that Eq. (15.31) is only approximately valid since there is a small contribution
from the kinetic energy term, which is however fully accounted for in an actual finite
temperature DFT calculation.
While the inclusion of F
el
0
is important for metals at realistic temperatures, for
semiconductors it is negligible except for narrow band gap semiconductors, where
also partial occupations f
i
may occur. In metals,
~
F
0
el
can become significant partic-
ularly if the density of states shows a peak close to the Fermi energy as found for
instance for d-states in Pt, Pd, Rh, or Ir [19].
15.2.2.3 Quasiharmonic Atomic Excitations
Neglecting for the moment the higher order terms in Eq. (15.29), the quasiharmonic
approximation results. The necessary steps to compute this contribution are
D
k;l
ðV; TÞ :¼
1
ffiffiffiffiffiffiffiffiffiffiffiffi
M
k
M
l
p
q
2
F
el
ðfR
I
g; V; TÞ
qR
k
qR
l
fR
0
I
g
! DðV; TÞw
i
¼ v
2
i
ðV; TÞw
i
! E
qh
fn
i
g
ðV; TÞ¼F
el
0
ðV; TÞþ
X
i
hv
i
ðV; TÞ
n
i
þ
1
2
;
ð15:32Þ
with n
i
the number of phonons in state i. In Eq. (15.32), the dynamical matrix D (with
elements D
k;l
) has been defined, which corresponds to the second derivative of F
el
scaled by the masses M
k
of the nuclei. The eigenvalue equation of D with eigenvectors
w
i
and eigenvalues v
2
i
defines the phonon frequencies v
i
. The energy E
qh
fn
i
g
for
a certain fixed phonon occupation configuration fn
i
g is given by a sum over the
frequencies weighted by the corresponding occupation numbers. Note the important
point that the phonon frequencies are not only volume dependent (quasiharmonic)
but also explicitly temperature dependent through the temperature dependence
of F
el
. This temperature dependence does not correspond to an anharmonic atomic
(i.e., phonon–phonon) interaction. To make this point explicit, we will use the
notation T
el
for this temperature, while the temperature determining the thermo-
dynamics of the nuclei will be denoted by T
nuc
. In fact, to speed up numeric
convergence, T
el
and T
nuc
can be varied independently of each other. Of course,
at the end of the calculation both have to be ensured to be equal to the actual external
temperature, i.e., T
el
¼ T
nuc
¼ T.
The final step of the quasiharmonic approximation is to approximate the eigen-
values
~
E
m
nuc
in Eq. (15.28) by the quasiharmonic energy E
qh
fn
i
g
which yields (the tilde in
15.2 Methodology
j
271
~
F indicating the approximation at this stage)
~
FðV; TÞ¼½
~
FðV; T
el
; T
nuc
Þ
T
el
¼T
nuc
¼T
; ð15:33Þ
with
~
F
ðV; T
el
; T
nuc
Þk
B
T ln
X
fn
i
g
e
b
~
E
qh
f
n
i
g
¼ … ¼ F
el
0
ðV; T
el
ÞþF
qh
ðV; T
el
; T
nuc
Þ;
ð15:34Þ
where the dots denote a series of straightforward transformations (see e.g., Ref. [26])
and with
F
qh
V; T
el
; T
nuc
¼
1
N
X
i
½
h
2
v
i
V; T
el
zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{
T
nuc
¼0K
zero point contr:
þ T
nuc
k
B
ln 1exp bhv
i
V; T
el
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
negativeðÞentropic contribution :
:
ð15:35Þ
Here, N is an appropriately chosen scaling factor which is, e.g., the number of
sampled frequencies divided by three (number of spatial dimensions) if F
qh
should
refer to a per atom quantity. Note that F
qh
fully contains quantum mechanical
effects (e.g., zero point vibrations) in contrast to a free energy obtained from a classical
molecular dynamics run.
The quasiharmonic equations presented above are general and equally well
applicable to defect and bulk cells. In practice, the calculation of the quasiharmonic
free energy contribution of the defect cell requires special attention to ensure
a consistent treatment with the corresponding perfect bulk cell. The reason for this
is the break of translational symmetry introduced by the presence of the defect. As
a consequence, the commonly applied [19] Fourier interpolation to generate dense
wave vector meshes in the Brillouin zone cannot be employed to the defect cell.
To nonetheless profit from the advantage of a Fourier interpolation for the perfect
bulk (e.g., a significantly improved description of the low temperature free energy)
a correction scheme as e.g., proposed in Ref. [18] should be used.
15.2.2.4 Anharmonic Atomic Excitations: Thermodynamic Integration
Let us now consider the higher order terms in Eq. (15.30). An intuitive and
straightforward approach to include these terms would be an ab initio based classical
molecular dynamics run. The fact that a direct computation of the free energy
employing conventional molecular dynamics is not feasible (the free energy is an
entropic quantity [27]) could be circumvented by calculating the inner energy and
integrating it with respect to temperature. It turns out that such a naive phase space
sampling leads to infeasibly long computational times [18]. Therefore, the develop-
ment and application of highly efficient sampling strategies to perform the ther-
modynamic averages is crucial.
272
j
15 Formation Energies of Point Defects at Finite Temperatures
A fundamental concept is thermodynamic integration. The key idea is to start from
a reference with an analytically known free energy or with a free energy that can be
obtained numerically extremely fast. This reference is coupled adiabatically to the
true ab initio potential energy surface and only the difference in free energies is
sampled. Such an approach will be efficient, if one manages to construct a reference
which closely approximates the true ab initio potential energy surface thus yielding
a small free energy difference.
A good reference to perform thermodynamic integration is the quasiharmonic
potential energy surface. Therefore, a system with free energy F
l
is introduced which
couples the electronic þ quasiharmonic system [having the free energy
~
F as given in
Eq. (15.34)] to the full system (having free energy F) by the adiabatic switching
parameter l. The boundary conditions are defined to be
½F
l
l¼1
¼ F and ½F
l
l ¼ 0
¼
~
F; ð15:36Þ
so that the anharmonic free energy is given by
F
ah
¼ F
~
F ¼½F
l
l¼1
½F
l
l¼0
: ð15:37Þ
In principle, methods based on the quantum-classical isomorphism exist allowing
to perform thermodynamic integration including quantum mechanical effects [28].
Unfortunately, their actual application requires very large computational efforts
and at present only investigations employing empirical potential energy surfaces
can be afforded [28]. In practice, a quantum mechanical treatment of F
ah
is not
expected to be necessary since the quasiharmonic free energy contains the major
part of quantum mechanical effects (as derived in the previous section). The
reason for this is that at low temperatures, where quantum effects are important,
the quasiharmonic approximation is excellent thus fully accounting for such
effects, while at high temperatures where the quasiharmonic approximation fails
and anharmonicity becomes significant quantum effects become small/negligible.
The anharmonic contribution can therefore be calculated in the classical limit
(superscript clas; note that correspondingly a classical quasiharmonic reference
is used):
F
clas;ah
:¼ F
clas
ðF
el
0
þF
clas;qh
Þ¼½F
clas
l
l¼1
½F
clas
l
l¼0
¼
ð
1
0
dl
qF
clas
l
ql
¼
#
Eq: ð15:39Þ
ð
1
0
dl
*
qF
el
l
ql
+
T;l
¼
#
ergodicity
hypothesis
ð
1
0
dl
*
qF
el
l
ql
+
t;l
:
ð15:38Þ
Here, F
el
l
is the l dependent electronic free energy surface determining the classical
motion of the nuclei in the coupled system, hi
T;l
denotes the thermodynamic, and
hi
t;l
the time average at a given l. To obtain the first equality in the third line, we have
used:
15.2 Methodology
j
273
qF
clas
l
ql
¼
1
Z
clas
l
ð
dR
I
V
ph
e
bF
el
l
ðfR
I
gÞ
qF
el
l
ql
¼
*
qF
el
l
ql
+
T;l
;
ð15:39Þ
with the standard definitions of a classical free energy F
clas
, partition function Z
clas
,
and phase space volume V
ph
. Besides the fixed boundary conditions, Eq. (15.36), any
type of coupled system can be chosen. In practice, a simple linear coupling to the
quasiharmonic reference (restoring the explicit notation for the volume and tem-
perature dependencies),
F
el
l
ðfR
I
g; V; T
el
Þ¼lF
el
ðfR
I
g; V; T
el
Þ
þð1 lÞ
F
el
0
ðV; T
el
Þþ
X
k;l
ffiffiffiffiffiffiffiffiffiffiffiffi
M
k
M
l
p
2
u
k
u
l
D
k;l
ðV; T
el
Þ
;
ð15:40Þ
yields computationally efficient results. Finally, the anharmonic free energy reads:
F
clas;ah
ðV; T
el
; T
nuc
Þ¼
ð
1
0
dl
F
el
ðV; T
el
Þ F
el
0
ðV; T
el
Þ
X
k;l
ffiffiffiffiffiffiffiffiffiffiffiffi
M
k
M
l
p
2
u
k
u
l
D
k;l
ðV; T
el
Þ
t;l
: ð15:41Þ
Note that the dependence of F
clas;ah
on T
nuc
is hidden in the time average
:
t
which
needs to be obtained e.g. from a molecular dynamics simulation. Combining the
various contributions, the free energy of the system is given by:
F ¼ E
el
g;0
þ
e
F
el
0
þF
qh
þF
clas;ah
¼½E
el
g;0
ðVÞþ
e
F
el
0
ðV; T
el
ÞþF
qh
ðV; T
el
; T
nuc
Þ
þF
clas;ah
ðV; T
el
; T
nuc
Þ
T
el
¼T
nuc
¼T
: ð15:42Þ
We stress that the dominant part of quantum mechanical effects contributing to F is
accounted for by F
qh
. For the case of (wide band gap) semiconductors the
influence of T
el
will be small (cf. Section 15.2.2.2) and the free energy can be
approximated by:
FðV; TÞ½E
el
g;0
ðVÞþF
qh
ðV; 0K; T
nuc
ÞþF
clas;ah
ðV; 0K; T
nuc
Þ
T
nuc
¼T
:
ð15:43Þ
15.2.2.5 Anharmonic Atomic Excitations: Beyond the Thermodynamic Integration
While the thermodynamic integration approach boosts the efficiency by a few orders
of magnitude as compared to conventional molecular dynamics simulations, the
calculation of a full free energy surface FðV; TÞfrom ab initio is still a formidable task
274
j
15 Formation Energies of Point Defects at Finite Temperatures
even on todays high performance CPU architecture. Therefore, there are currently
active efforts at exploring new methods to reduce computational time [18, 29, 30].
Wu and Wentzcovitch [29], e.g., developed a semi empirical ansatz based on
a single adjustable parameter to describe anharmonicity. To explain their approach, it
is useful to note first that already the quasiharmonic frequencies v
qh
i
contain an
implicit temperature dependence due to thermal expansion: v
qh
i
¼ v
qh
i
ðV
eq
ðTÞÞ.
Here, V
eq
ðTÞ denotes the thermal expansion, i.e., the equilibrium volume as
a function of temperature T (for a fixed pressure). The thermal expansion is itself
determined by the v
qh
i
so the problem must be solved self consistently.
The ansatz of Wu and Wentzcovitch [29] is to modify the implicit temperature
dependence. Specifically, the following transformation to a renormalized frequency
v
renorm
i
is proposed:
v
qh
i
V
eq
ðTÞ
!v
renorm
i
ðTÞ :¼ v
qh
i
V
0
ðTÞ
; ð15:44Þ
with a modified thermal expansion
V
0
ðTÞ¼V
eq
ðTÞ 1C
V
eq
ðTÞV
eq
0
V
eq
0
; ð15:45Þ
where C is an adjustable parameter and V
eq
0
¼ V
eq
ðT ¼ 0KÞ. The renormalized
frequencies v
renorm
i
replace then the quasiharmonic frequencies in the quasihar-
monic free energy expression Eq. (15.35). The resulting free energy surface includes
therefore anharmonicity in an approximative manner and can be used to derive
anharmonic thermodynamic quantities such as thermal expansion or defect con-
centrations. In the original work [29], C is adjusted by fitting to experimental data.
In a more recent work [30], it was shown how the parameter can be obtained using
ab initio based thermodynamic integration. The main advantage of the approach is
its effectiveness, since to determine C and thus to generate a complete anharmonic
free energy surface only a single thermodynamic integration run at fixed V and T
suffices. The disadvantage is however that the v
renorm
i
(and thus the anharmonic
free energy) are assumed to have a certain temperature dependence, specifically that
they scale linearly with the quasiharmonic thermal expansion. This assumption has
been shown to work well for diamond and MgO [30] up to high temperatures, it
cannot be however a priori assumed in general and further evaluations are
necessary.
An alternative approach that is able to account for the full anharmonic free energy
surface without assuming any specific temperature dependence was developed
recently in Ref. [18] and termed upsampled thermodynamic integration using
Langevin dynamics (UP-TILD) method. The key ideas are as follows. In a first step,
a usual thermodynamic integration run for discrete values of l is performed:
*
qF
el
l
ql
+
low
t;l
¼
*
F
el
ðV; T
el
ÞF
el
0
ðV; T
el
Þ
X
k;l
ffiffiffiffiffiffiffiffiffiffiffiffi
M
k
M
l
p
2
u
k
u
l
D
k;l
ðV; T
el
Þ
+
low
t;l
:
ð15:46Þ
15.2 Methodology
j
275
A critical insight gained in Ref. [18] is that instead of employing fully converged DFT
parameters (e.g., with respect to basis set or electronic k sampling) computationally
much less demanding DFT parameters can be used (indicated by the superscript
low). In particular, they need to be chosen such that the corresponding phase space
distribution (termed R
I
fg
low
t
in the following) closely resembles the phase space
distribution R
I
fg
high
t
as would be obtained from fully/highly converged parameters.
As a consequence of using the low converged parameters, the thermodynamic
integration is computationally very ef ficient and speed up factors of 30 are
achievable. Applying this approach makes it therefore possible to sample various
l, V, and T values even on modest computer resources. The resulting free energy
surface, i.e., hqF
el
l
/qli
low
t;l
, needs however to be corrected in a second step. The actual
correction step is straightforward, provided that the above stated condition R
I
fg
low
t
close to R
I
fg
high
t
holds. In such a case, a typically small set of N
UP
uncorrelated
structures R
I
fg
low
t
(indexed with t
u
) is extracted from R
I
fg
low
t
and the upsampling
average hDF
el
i
UP
l
is calculated as:
DF
el
hi
UP
l
¼
1
N
UP
X
N
UP
u
F
el;low
R
I
fg
low
t
u
F
el;low
R
0
I
F
el;high
R
I
fg
low
t
u
F
el;high
R
0
I
hi
:
ð15:47Þ
Here,F
el;low
(F
el;high
) refers to the electronic free energy calculated using the low (high)
converged set of DFT parameters. The l dependence of hD F
el
i
UP
l
is hidden in the
trajectory R
I
fg
low
t
, which is additionally dependent on the volume and temperature. In
the last step, the quantity of interest, the converged hqF
el
l
/qli
high
t;l
is obtained
qF
el
l
/ql
high
t;l
¼ qF
el
l
/ql
high
t;l
DF
el
UP
l
;
and thus the anharmonic free energy:
F
clas;ah
¼
ð
1
0
dlhqF
el
l
/qli
high
t;l
: ð15:48Þ
The efficiency of the UP-TILD method is a direct consequence of the fact that the
upsamling average hDF
el
i
UP
l
, which involves the computation of the CPU time con-
sumingwellconvergedDFTfreeenergiesF
el;high
,convergesextremelyfastwithrespectto
the number of uncorrelated configurations. In practice, less than 100 configurations are
found to be sufficient to achieve an accuracy of better than 1 meV/atom in the free
energy [18]. This small number of configurations has to be compared with the number of
moleculardynamics steps needed to convergea thermodynamic integrationrun whichis
in the rangeof several thousands, i.e., two orders of magnitude larger. Another advantage
is that in many cases the l dependence is invariant with respect to the upsampling
procedure allowing to calculate the shift in Eq. (15.48) for a single l value only. Using this
invariance reduces computational costs further as illustrated in Figure 15.4.
The various methods needed to compute the electronic, quasiharmonic, and
anharmonic contributions to the free energy surface are summarized in Figure 15.5.
276
j
15 Formation Energies of Point Defects at Finite Temperatures
Figure 15.4 (online colour at: www.pss-b.com) CPU time reductions in calculating the anharmonic free energy as reported for the example of bulk
aluminum including point defects (32 atomic fcc cell) [18]; MD: molecular dynamics; TI: thermodynamic integration; UP-TILD: upsampled TI using Langevin dynamics
(see text for details).
15.2 Methodology
j
277