of 17
Perturbo
: a software package for
ab initio
electron-phonon interactions,
charge transport and ultrafast dynamics
Jin-Jian Zhou
a,
, Jinsoo Park
a
, I-Te Lu
a
, Ivan Maliyov
a
, Xiao Tong
a
, Marco Bernardi
a,
a
Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA.
Abstract
Perturbo
is a software package for first-principles calculations of charge transport and ultrafast carrier dynamics
in materials. The current version focuses on electron-phonon interactions and can compute phonon-limited transport
properties such as the conductivity, carrier mobility and Seebeck coefficient. It can also simulate the ultrafast nonequi-
librium electron dynamics in the presence of electron-phonon scattering.
Perturbo
uses results from density functional
theory and density functional perturbation theory calculations as input, and employs Wannier interpolation to reduce
the computational cost. It supports norm-conserving and ultrasoft pseudopotentials, spin-orbit coupling, and polar
electron-phonon corrections for bulk and 2D materials. Hybrid MPI plus OpenMP parallelization is implemented to
enable efficient calculations on large systems (up to at least 50 atoms) using high-performance computing. Taken to-
gether,
Perturbo
provides efficient and broadly applicable
ab initio
tools to investigate electron-phonon interactions
and carrier dynamics quantitatively in metals, semiconductors, insulators, and 2D materials.
Keywords:
Charge transport, Ultrafast dynamics, Electron-phonon interactions, Wannier functions
PROGRAM SUMMARY
Program Title:
Perturbo
Program summary URL:
https://perturbo-code.github.io
Licensing provisions:
GNU General Public Licence 3.0
Programming language:
Fortran
,
Python
External routines/libraries:
LAPACK
,
HDF5
,
MPI
,
OpenMP
,
FFTW,
Quantum-ESPRESSO
Nature of problem:
Computing transport properties from first-
principles in materials, including the electrical conductivity,
carrier mobility and Seebeck coefficient; Simulating ultrafast
nonequilibrium electron dynamics, such as the relaxation of
excited carriers via interactions with phonons.
Solution method:
We implement the first-principles Boltzmann
transport equation, which employs materials properties such as
the electronic structure, lattice dynamics, and electron-phonon
collision terms computed with density functional theory and
density functional perturbation theory. The Boltzmann trans-
port equation is solved numerically to compute charge trans-
port and simulate ultrafast carrier dynamics. Wannier interpo-
lation is employed to reduce the computational cost.
Additional comments:
Hybrid MPI plus OpenMP paralleliza-
tion is implemented to run large calculations and take advan-
tage of high-performance computing. Most results are output
to HDF5 file format, which is portable and convenient for post-
processing using high-level languages such as Python and Ju-
lia.
Corresponding author.
E-mail addresses:
jjzhou@caltech.edu (J.-J. Zhou), bmarco@caltech.edu (M. Bernardi).
1. Introduction
Understanding the dynamical processes involving elec-
trons, lattice vibrations (phonons), atomic defects, and
photons in the solid state is key to developing the next
generation of materials and devices [1, 2]. Due to the in-
creasing complexity of functional materials, there is a criti-
cal need for computational tools that can take into account
the atomic and electronic structure of materials and make
quantitative predictions on their physical properties. The
vision behind
Perturbo
is to provide a unified platform
and a validated code that can be applied broadly to com-
pute the interactions, transport and ultrafast dynamics
of electrons and excited states in materials [3]. The goal
is to facilitate basic scientific discoveries in materials and
devices by advancing microscopic understanding of carrier
dynamics, while creating a sustainable software element
able to address the demands of the computational physics
community.
Perturbo
builds on established first-principles meth-
ods. It uses density functional theory (DFT) and den-
sity functional perturbation theory (DFPT) [4] as a start-
ing point for computing electron dynamics. It reads the
output of DFT and DFPT calculations, for now from the
Quantum Espresso
(QE) code [5, 6], and uses this data
to compute electron interactions, charge transport and
ultrafast dynamics. The current distribution focuses on
electron-phonon (
e
-ph) interactions and the related phonon-
limited transport properties [7], including the electrical
conductivity, mobility and the Seebeck coefficient. It can
also simulate the ultrafast nonequilibrium electron dynam-
Preprint submitted to Computer Physics Communications
February 7, 2020
arXiv:2002.02045v1 [cond-mat.mtrl-sci] 6 Feb 2020
ics in the presence of
e
-ph interactions. The developer
branch, which is not publicly available yet, also includes
routines for computing spin [8], electron-defect [9, 10],
and electron-photon interactions [11], as well as advanced
methods to compute the ultrafast dynamics of electrons
and phonons in the presence of electric and magnetic fields.
These additional features will be made available in future
releases.
The transport module of
Perturbo
enables accurate
calculations of charge transport in a wide range of func-
tional materials. In its most basic workflow,
Perturbo
computes the conductivity and mobility as a function of
temperature and carrier concentration, either within the
relaxation time approximation (RTA) or with an iterative
solution of the linearized Boltzmann transport equation
(BTE) [12, 13]. The ultrafast dynamics module explic-
itly evolves in time the electron BTE (while keeping the
phonon occupations fixed), enabling investigations of the
ultrafast electron dynamics starting from a given initial
electron distribution [14]. Our routines can carry out these
calculations in metals, semiconductors, insulators, and 2D
materials. An efficient implementation of long-range
e
-ph
interactions is employed for polar bulk and 2D materials.
Materials with spin-orbit coupling (SOC) are treated us-
ing fully relativistic pseudopotentials [8, 13]. Both norm-
conserving and ultrasoft pseudopotentials are supported.
Quantities related to
e
-ph interactions can be easily ob-
tained, stored and analyzed.
Perturbo
is implemented in modern Fortran with a
modular code design. All calculations employ intuitive
workflows. The code is highly efficient thanks to its hy-
brid MPI (Message Passing Interface) and OpenMP (Open
Multi-Processing) parallelization. It can run on record-
large unit cells with up to at least 50 atoms [15], and
its performance scales up to thousands of CPU cores. It
conveniently writes files using the HDF5 format, and is
suitable for both high-performance supercomputers and
smaller computer clusters.
Target users include both experts in first-principles cal-
culations and materials theory as well as experimental re-
searchers and teams in academic or national laboratories
investigating charge transport, ultrafast spectroscopy, ad-
vanced functional materials, and semiconductor or solid-
state devices.
Perturbo
will equip these users with an
efficient quantitative tool to investigate electron interac-
tions and dynamics in broad families of materials, filling a
major void in the current software ecosystem.
The paper is organized as follows: Sec. 2 discusses the
theory and numerical methods implemented in the code;
Sec. 3 describes the code capabilities and workflows; Sec. 4
delves deeper into selected technical aspects; Sec. 5 shows
several example calculations provided as tutorials in the
code; Sec. 6 discusses the parallelization strategy and the
scaling of the code on high performance supercomputers.
We conclude in Sec. 7 by summarizing the main features
and planned future development of
Perturbo
.
2. Methodology
2.1. Boltzmann transport equation
The current release of
Perturbo
can compute charge
transport and ultrafast dynamics in the framework of the
semiclassical BTE. The BTE describes the flow of the elec-
tron occupations
f
n
k
(
r
,t
) in the phase-space variables of
relevance in a periodic system, the crystal momentum
k
and spatial coordinate
r
:
∂f
n
k
(
r
,t
)
∂t
=
[
r
f
n
k
(
r
,t
)
·
v
n
k
+
~
1
k
f
n
k
(
r
,t
)
·
F
]
+
I
[
f
n
k
]
,
(1)
where
n
is the band index and
v
n
k
are band velocities. The
time evolution of the electron occupations is governed by
the so-called drift term due to external fields
F
and the
collision term
I
[
f
n
k
], which captures electron scattering
processes due to phonons or other mechanisms [16]. In
Perturbo
, the fields are assumed to be slowly varying
and the material homogeneous, so
f
n
k
does not depend on
the spatial coordinates and its spatial dependence is not
computed explicitly.
The collision integral
I
[
f
n
k
] is a sum over a large num-
ber of scattering processes in momentum space, and it is
very computationally expensive because it involves Bril-
louin zone (BZ) integrals on fine grids. Most analytical and
computational treatments simplify the scattering integral
with various approximations. A common one is the RTA,
which assumes that the scattering integral is proportional
to the deviation
δf
n
k
of the electron occupations from the
equilibrium Fermi-Dirac distribution,
I
[
f
n
k
] =
δf
n
k
;
the relaxation time
τ
is either treated as a constant empir-
ical parameter [17, 18] or as a state-dependent quantity,
τ
n
k
.
Perturbo
implements the first-principles formalism
of the BTE, which employs materials properties obtained
with quantum mechanical approaches, using the atomic
structure of the material as the only input. The elec-
tronic structure is computed using DFT and the lattice
dynamical properties using DFPT. The scattering integral
is computed only for
e
-ph processes in the current release,
while other scattering mechanisms, such as electron-defect
and electron-electron scattering, are left for future releases.
The scattering integral due to
e
-ph processes can be writ-
ten as
I
e
ph
[
f
n
k
] =
2
π
~
1
N
q
m
q
ν
|
g
mnν
(
k
,
q
)
|
2
×
[
δ
(
ε
n
k
~
ω
ν
q
ε
m
k
+
q
)
×
F
em
+
δ
(
ε
n
k
+
~
ω
ν
q
ε
m
k
+
q
)
×
F
abs
]
,
(2)
where
N
q
is the number of
q
-points used in the summation,
and
g
mnν
(
k
,
q
) are
e
-ph matrix elements (see Sec. 2.2.3)
quantifying the probability amplitude for an electron to
scatter from an initial state
|
n
k
to a final state
|
m
k
+
q
,
by emitting or absorbing a phonon with wavevector
q
and
2
mode index
ν
; here and below,
ε
n
k
and
~
ω
ν
q
are the en-
ergies of electron quasiparticles and phonons, respectively.
The phonon absorption (
F
abs
) and emission (
F
em
) terms
are defined as
F
abs
=
f
n
k
(1
f
m
k
+
q
)
N
ν
q
f
m
k
+
q
(1
f
n
k
) (
N
ν
q
+ 1)
,
F
em
=
f
n
k
(1
f
m
k
+
q
) (
N
ν
q
+ 1)
f
m
k
+
q
(1
f
n
k
)
N
ν
q
(3)
where
N
ν
q
are phonon occupations.
2.1.1. Ultrafast carrier dynamics
Perturbo
can solve the BTE numerically to simulate
the evolution in time
t
of the electron occupations
f
n
k
(
t
)
due to
e
-ph scattering processes, starting from an initial
nonequilibrium distribution. In the presence of a slowly
varying external electric field, and assuming the material
is homogeneous, we rewrite the BTE in Eq. (1) as
∂f
n
k
(
t
)
∂t
=
e
E
~
·∇
k
f
n
k
(
t
) +
I
e
ph
[
f
n
k
]
,
(4)
where
e
is the electron charge and
E
the external electric
field. We solve this non-linear integro-differential equation
numerically (in the current release, only for
E
= 0), us-
ing explicit time-stepping with the 4th-order Runge-Kutta
(RK4) or Euler methods. The RK4 solver is used by de-
fault due to its superior accuracy. The Euler method is
much faster than RK4, but it is only first-order accurate
in time, so it should be tested carefully and compared
against RK4.
Starting from an initial nonequilibrium electron distri-
bution
f
n
k
(
t
0
), we evolve in time Eq. (4) with a small time
step (typically in the order of 1 fs) to obtain
f
n
k
(
t
) as a
function of time. One application of this framework is to
simulate the equilibration of excited electrons [14, 19], in
which we set
E
= 0 and simulate the time evolution of the
excited electron distribution as it approaches its equilib-
rium Fermi-Dirac value through
e
-ph scattering processes.
The phonon occupations are kept fixed (usually, to their
equilibrium value at a given temperature) in the current
version of the code. Another application, which is cur-
rently under development, is charge transport in high elec-
tric fields.
2.1.2. Charge transport
In an external electric field, the drift and collision terms
in the BTE balance out at long enough times
the field
drives the electron distribution out of equilibrium, while
the collisions tend to restore equilibrium. At steady state,
a nonequilibrium electron distribution is reached, for which
∂f
n
k
/∂t
= 0. The BTE for transport at steady state be-
comes
e
E
~
·∇
k
f
n
k
(
t
) =
I
e
ph
[
f
n
k
]
.
(5)
When the electric field is relatively weak, the steady-state
electron distribution deviates only slightly from its equi-
librium value. As is usual, we expand
f
n
k
around the
equilibrium Fermi-Dirac distribution,
f
0
n
k
, and keep only
terms linear in the electric field:
f
n
k
=
f
0
n
k
+
f
1
n
k
+
O
(
E
2
)
=
f
0
n
k
+
e
E
·
F
n
k
∂f
0
n
k
∂ε
n
k
+
O
(
E
2
)
,
(6)
where
F
n
k
characterizes the first-order deviation from equi-
librium of the electron distribution. We substitute Eq. (6)
into both sides of Eq. (5), and obtain a linearized BTE for
the distribution deviation
F
n
k
keeping only terms up to
first-order in the electric field:
F
n
k
=
τ
n
k
v
n
k
+
τ
n
k
N
q
m,ν
q
F
m
k
+
q
W
ν
q
n
k
,m
k
+
q
,
(7)
where
τ
n
k
is the electron relaxation time, computed as the
inverse of the scattering rate,
τ
n
k
= Γ
1
n
k
. The scattering
rate Γ
n
k
is given by
Γ
n
k
=
1
N
q
m,ν
q
W
ν
q
n
k
,m
k
+
q
.
(8)
The scattering probability
W
ν
q
n
k
,m
k
+
q
involves both phonon
emission and absorption processes:
W
ν
q
n
k
,m
k
+
q
=
2
π
~
|
g
mnν
(
k
,
q
)
|
2
×
[
δ
(
ε
n
k
~
ω
ν
q
ε
m
k
+
q
)
(
1 +
N
0
ν
q
f
0
m
k
+
q
)
+
δ
(
ε
n
k
+
~
ω
ν
q
ε
m
k
+
q
)
(
N
0
ν
q
+
f
0
m
k
+
q
)
]
.
(9)
where
N
0
ν
q
are equilibrium Bose-Einstein phonon occupa-
tions. Note that since
τ
n
k
is an electron quasiparticle life-
time, it can be written equivalently as the imaginary part
of the
e
-ph self-energy [20],
τ
1
n
k
= 2ImΣ
e
ph
n
k
/
~
.
In the RTA, we neglect the second term in Eq. (7) and
obtain
F
nk
=
τ
n
k
v
n
k
. In some cases, the second term
in Eq. (7) cannot be neglected. In metals, a commonly
used scheme to approximate this term is to add a factor
of (1
cos
θ
k
,
k
+
q
) to Eq. (9), where
θ
k
,
k
+
q
is the scatter-
ing angle between
k
and
k
+
q
. The resulting so-called
“transport relaxation time” is then used to compute the
transport properties [20].
Perturbo
implements a more
rigorous approach and directly solves Eq. (7) using an it-
erative method [21, 22], for which we rewrite Eq. (7) as
F
i
+1
n
k
=
F
0
nk
+
τ
n
k
N
q
m,ν
q
F
i
m
k
+
q
W
ν
q
n
k
,m
k
+
q
.
(10)
In the iterative algorithm, we choose in the first step
F
0
nk
=
τ
n
k
v
n
k
, and then compute the following steps using Eq. (10)
until the difference
|
F
i
+1
n
k
F
i
n
k
|
is within the convergence
threshold.
Once
F
nk
has been computed, either within the RTA
or with the iterative solution of the BTE in Eq. (10), the
conductivity tensor is obtained as
σ
αβ
=
1
E
β
·
S
N
k
n
k
e
v
α
n
k
·
f
1
n
k
=
e
2
S
N
k
n
k
v
α
n
k
F
β
n
k
(
∂f
0
n
k
∂ε
n
k
)
,
(11)
3
where
α
and
β
are Cartesian directions, Ω is the volume
of the unit cell, and
S
is the spin degeneracy. We also
compute the carrier mobility tensor,
μ
αβ
=
σ
αβ
/
(
en
c
), by
dividing the conductivity tensor through the carrier con-
centration
n
c
.
In our implementation, we conveniently rewrite Eq. (11)
as
σ
αβ
=
e
2
dE
(
∂f
0
/∂E
αβ
(
E
)
,
(12)
where Σ
αβ
(
E
) is the transport distribution function (TDF)
at energy
E
,
Σ
αβ
(
E
) =
S
N
k
n
k
v
α
n
k
F
β
nk
δ
(
E
ε
n
k
)
,
(13)
which is computed in
Perturbo
using the tetrahedron
integration method [23]. The integrand in Eq. (12) can
be used to characterize the contributions to transport as
a function of electron energy [12, 13]. The code can also
compute the Seebeck coefficient
S
from the TDF, using
[
σS
]
αβ
=
e
T
dE
(
∂f
0
/∂E
)(
E
μ
αβ
(
E
)
,
(14)
where
μ
is the chemical potential and
T
is the temperature.
2.2. Electrons, phonons, and
e
-ph interactions
As discussed above, the electronic structure and phonon
dispersion are computed with DFT and DFPT. In princi-
ple, the
e
-ph matrix elements can also be computed with
these methods and used directly for carrier dynamics cal-
culations. However, to converge transport and ultrafast
dynamics with the BTE, the
e
-ph matrix elements and the
scattering processes need to be computed on ultra-dense
k
-
and
q
-point BZ grids with roughly 100
×
100
×
100 or more
points. Therefore, the computational cost is prohibitive for
direct DFPT calculations, and we resort to interpolation
techniques to obtain the
e
-ph matrix elements and other
relevant quantities on fine grids, starting from DFT and
DFPT calculations on coarser BZ grids, typically of order
10
×
10
×
10.
2.2.1. Wannier interpolation of the electronic structure
We use the Wannier interpolation to compute efficiently
the electron energy and band velocity on ultra-fine
k
-point
grids [24]. We first perform DFT calculations on a regular
coarse grid with points
k
c
, and obtain the electron energies
ε
n
k
c
and Bloch wavefunctions
|
ψ
n
k
c
. We construct maxi-
mally localized Wannier functions
|
n
R
e
, with index
n
and
centered in the cell at
R
e
, from the Bloch wavefunctions
using the
Wannier90
code [25, 26]:
|
n
R
e
=
1
N
e
m
k
c
e
i
k
c
·
R
e
U
mn
(
k
c
)
|
ψ
m
k
c
,
(15)
where
N
e
is the number of
k
c
-points in the coarse grid,
and
U
(
k
c
) are the unitary matrices transforming the Bloch
wavefunctions to a Wannier gauge [27],
|
ψ
(
W
)
n
k
c
=
m
U
mn
(
k
c
)
|
ψ
m
k
c
.
(16)
For entangled band structures,
U
(
k
c
) are not in general
square matrices since they are also used to extract a smooth
subspace from the original DFT Bloch eigenstates [28].
We compute the electron Hamiltonian in the Wannier
function basis,
H
nn
(
R
e
) =
n
0
|
ˆ
H
|
n
R
e
=
1
N
e
k
c
e
i
k
c
·
R
e
[
U
(
k
c
)
H
(
k
c
)
U
(
k
c
)
]
nn
(17)
where
H
(
k
c
) is the Hamiltonian in the DFT Bloch eigen-
state basis,
H
nm
(
k
c
) =
ε
n
k
c
δ
nm
. The Hamiltonian in the
Wannier basis,
H
nn
(
R
e
), can be seen as an
ab initio
tight-
binding model, with hopping integrals from the Wannier
orbital
n
in the cell at
R
e
to the Wannier orbital
n
in the
cell at the origin. Due to the localization of the Wannier
orbitals, the hopping integrals decay rapidly with
|
R
e
|
, so
a small set of
R
e
vectors is sufficient to represent the elec-
tronic structure of the system.
Starting from
H
nn
(
R
e
), we obtain the band energy
ε
n
k
and band velocity
v
n
k
at any desired
k
-point. We first
compute the Hamiltonian matrix
H
(
W
)
(
k
) in the basis of
Bloch sums of Wannier functions using an inverse discrete
Fourier transform, and then diagonalize it through a uni-
tary rotation matrix
U
(
k
) satisfying
H
(
W
)
(
k
) =
R
e
e
i
k
·
R
e
H
(
R
e
) =
U
(
k
)
H
(
H
)
(
k
)
U
(
k
)
,
(18)
where
H
(
H
)
nm
(
k
) =
ε
n
k
δ
nm
, and
ε
n
k
and
U
(
k
) are the eigen-
values and eigenvectors of
H
(
W
)
(
k
), respectively. One can
also obtain the corresponding interpolated Bloch eigen-
states as
|
ψ
n
k
=
m
U
mn
(
k
)
|
ψ
(
W
)
m
k
=
m
U
mn
(
k
)
R
e
e
i
k
·
R
e
|
m
R
e
.
(19)
The band velocity in the Cartesian direction
α
is computed
as
~
v
α
n
k
= [
U
(
k
)
H
(
W
)
α
(
k
)
U
(
k
)]
nn
,
(20)
where
H
(
W
)
α
(
k
) is the
k
-derivative of
H
(
W
)
(
k
) in the
α
-
direction, evaluated analytically using
H
(
W
)
α
(
k
) =
α
H
(
W
)
(
k
) =
R
e
e
i
k
·
R
e
H
(
R
e
)
·
(
iR
α
e
)
.
(21)
An appropriate extension of Eq. (20) is used for degenerate
states [29].
2.2.2. Interpolation of the phonon dispersion
The lattice dynamical properties are first obtained us-
ing DFPT on a regular coarse
q
c
-point grid. Starting
4
from the dynamical matrices
D
(
q
c
), we compute the inter-
atomic force constants (IFCs)
D
(
R
p
) (here, without the
mass factor) through a Fourier transform [30, 31],
D
(
R
p
) =
1
N
p
q
c
e
i
q
c
·
R
p
D
(
q
c
)
.
(22)
where
N
p
is the number of
q
c
-points in the coarse grid.
If the IFCs are short-ranged, a small set of
D
(
R
p
), ob-
tained from dynamical matrices on a coarse
q
c
-point grid,
is sufficient to obtain the dynamical matrix at any desired
q
-point with an inverse Fourier transform,
D
(
q
) =
R
p
e
i
q
·
R
p
D
(
R
p
)
.
(23)
We obtain the phonon frequencies
ω
ν
q
and displacement
eigenvectors
e
ν
q
by diagonalizing
D
(
q
).
2.2.3. Wannier interpolation of the
e
-ph matrix elements
The key quantities for
e
-ph scattering are the
e
-ph ma-
trix elements
g
mnν
(
k
,
q
), which appear above in Eq. (2).
They are given by
g
mnν
(
k
,
q
) =
~
2
ω
ν
q
κα
e
κα
ν
q
M
κ
ψ
m
k
+
q
|
q
,κα
V
|
ψ
n
k
,
(24)
where
|
ψ
n
k
and
|
ψ
m
k
+
q
are the wavefunctions of the ini-
tial and final Bloch states, respectively, and
q
,κα
V
is the
perturbation potential due to lattice vibrations, computed
as the variation of the Kohn-Sham potential
V
with re-
spect to the atomic displacement of atom
κ
(with mass
M
κ
) in the Cartesian direction
α
:
q
,κα
V
=
R
p
e
i
q
·
R
p
∂V
∂R
pκα
.
(25)
We obtain this perturbation potential as a byproduct of
the DFPT lattice dynamical calculations at a negligible
additional cost.
We compute the bra-ket in Eq. (24) directly, using the
DFT Bloch states on a coarse
k
c
-point grid and the per-
turbation potentials on a coarse
q
c
-point grid,
̃
g
κα
mn
(
k
c
,
q
c
) =
ψ
m
k
c
+
q
c
|
q
c
,κα
V
|
ψ
n
k
c
,
(26)
from which we obtain the
e
-ph matrix elements in the
Wannier basis [32, 33], ̃
g
κα
ij
(
R
e
,
R
p
), by combining Eq. (15)
and the inverse transformation of Eq. (25):
̃
g
κα
ij
(
R
e
,
R
p
) =
i
0
∂V
∂R
p,κα
j
R
e
=
1
N
e
N
p
k
c
,
q
c
e
i
(
k
c
·
R
e
+
q
c
·
R
p
)
̃
g
κα,
(
W
)
ij
(
k
c
,
q
c
)
,
(27)
where
̃
g
κα,
(
W
)
(
k
c
,
q
c
) =
U
(
k
c
+
q
c
)
̃
g
κα
(
k
c
,
q
c
)
U
(
k
c
)
are the matrix elements in the Wannier gauge. Similar to
the electron Hamiltonian in the Wannier basis, ̃
g
κα
ij
(
R
e
,
R
p
)
can be seen as a hopping integral between two localized
Wannier functions, one at the origin and one at
R
e
, due
to a perturbation caused by an atomic displacement at
R
p
.
If the interactions are short-ranged in real space, ̃
g
decays
rapidly with
|
R
e
|
and
|
R
p
|
, and computing it on a small
set of (
R
e
,
R
p
) lattice vectors is sufficient to fully describe
the coupling between electrons and lattice vibrations.
The
e
-ph matrix elements at any desired pair of
k
-
and
q
-points can be computed efficiently using the inverse
transformation in Eq. (27),
̃
g
κα
mn
(
k
,
q
) =
i,j
U
mi
(
k
+
q
)
U
jn
(
k
)
×
R
e
,
R
p
e
i
(
k
·
R
e
+
q
·
R
p
)
̃
g
κα
ij
(
R
e
,
R
p
)
,
(28)
where
U
(
k
) is the matrix used to interpolate the Bloch
states in Eq. (19).
The main requirement of this interpolation approach is
that the
e
-ph interactions are short-ranged and the
e
-ph
matrix elements in the local basis decay rapidly. There-
fore, the
e
-ph interpolation works equally well with local-
ized orbitals other than Wannier functions, as we have
shown recently using atomic orbitals [34]. The atomic or-
bital
e
-ph interpolation is still being extended and tested
in the
Perturbo
code, so the current release includes only
Wannier interpolation routines.
2.3. Polar corrections for phonons and
e
-ph interactions
The assumption that the IFCs and
e
-ph interactions
are short-ranged does not hold for polar semiconductors
and insulators. In polar materials, the displacement of ions
with a non-zero Born effective charge creates dynamical
dipoles, and the long-wavelength longitudinal optical (LO)
phonon mode induces a macroscopic electric field [30]. The
dipole-dipole interactions introduce long-range contribu-
tions to the IFCs and dynamical matrices [35], resulting in
the well-known LO-TO splitting in the phonon dispersion
at
q
0. For this reason, the dynamical matrix inter-
polation scheme in Eqs. (22)-(23) cannot provide correct
phonon dispersions at small
q
for polar materials. To ad-
dress this issue, a polar correction is typically used [31], in
which the dynamical matrix is separated into two contri-
butions, a short-range part that can be interpolated using
the Fourier transformation in Eqs. (22)-(23), and a long-
range part evaluated directly using an analytical formula
involving the Born effective charges and the dielectric ten-
sor [31].
Perturbo
implements this standard approach
to include the LO-TO splitting in polar materials.
The field due to the dynamical dipoles similarly intro-
duces long-range
e
-ph contributions
in particular, the
Fr ̈ohlich interaction [36], a long-range coupling between
electrons and LO phonons. The strength of the Fr ̈ohlich
e
-ph interaction diverges as 1
/q
for
q
0 in bulk materi-
als. As a result, the Wannier interpolation is impractical
5