Before running electron dynamics calculations using pertubo.x
, the user needs to carry out electronic and phonon calculations, with DFT and DFPT respectively. At present, Perturbo can read the output of DFT and DFPT calculations done with Quantum Espresso (QE). Once the relevant output files have been obtained from QE and Wannier90 (W90), the first step is to use qe2pert.x
to compute the eph matrix elements on a coarse \(\mathbf{k}\) and \(\mathbf{q}\) point Brillouin zone grid, to obtain eph matrix elements in Wannier function basis, and to store the data into the HDF5 format for perturbo.x
to read. The generation of this HDF5 file, called ‘prefix’_epr.h5, is discussed in this section of the manual.
The preparation stage consists of five steps:
 Run a selfconsistent (SCF) DFT calculation
 Run a phonon calculation using DFPT
 Run a nonSCF (nSCF) DFT calculation
 Run Wannier90 to obtain Wannier functions
 Run
qe2pert.x
In the following, we use silicon as an example. The input files for QE and W90 are in the directory “example02siliconqe2pert/pwphwann” of the repository of perturboexampleslight. As a reference, we also provide the results in a directory called “References”.
Step 1: SCF calculation
Here is the scf.in QE input file:
&CONTROL
prefix = 'si'
calculation = 'scf'
wf_collect = .true.
outdir='./tmp'
pseudo_dir='../pseudo'
/
&SYSTEM
ibrav = 2
celldm(1) = 10.264
nat = 2
ntyp = 1
ecutwfc = 40.0
/
&ELECTRONS
conv_thr = 1.0d15
mixing_mode = 'plain'
mixing_beta = 0.7
diagonalization = 'david'
diago_full_acc = .true.
/
&IONS
/
&CELL
press_conv_thr = 0.01
/
ATOMIC_SPECIES
Si 28.085 Si_DOJO_PBEsol.upf
ATOMIC_POSITIONS crystal
Si 0.00000000 0.00000000 0.00000000
Si 0.25000000 0.75000000 0.25000000
K_POINTS (automatic)
20 20 20 0 0 0
Run an SCF calculation (pw.x
QE executable) and obtain the QE ‘prefix’.save directory. In this case, we obtain ./tmp/si.save, which is needed for phonon and nSCF calculations.
Step 2: phonon calculation
Here is the ph.in input file:
Phonons on a uniform grid
&inputph
verbosity='debug'
tr2_ph=1.0d17
prefix='si'
ldisp=.true.
epsil=.true.
lqdir=.true.
outdir='./tmp'
fildyn = 'si.dyn.xml'
fildvscf = 'dvscf'
nq1=8, nq2=8, nq3=8,
/
Copy the QE ‘prefix’.save directory from the SCF run to the current directory:
$ cp r ../scf/tmp ./
In this step, make sure that the number of \(\mathbf{q}\) points is commensurate with the number of \(\mathbf{k}\) points used in the nSCF and Wannierization calculations. For example, a \(\mathbf{q}\) grid of 8x8x8 can be used with a wannierization \(\mathbf{k}\) grid of 8x8x8 or 16x16x16, but not with a 10x10x10 or 12x12x12 grid.
start_q
and last_q
parameters is also possible.Run the PHonon QE calculation (ph.x
QE executable).
Once the calculations are done, we collect the phonon data into a directory called save, created by running the shell script phcollect.sh.
$ ./phcollect.sh
or
$ bash phcollect.sh
The save directory contains all the information needed for PERTURBO to interface with QE. It includes the dynamical matrix files, phonon perturbation potentials, and the patterns.
Since phonon calculations using DFPT can be computationally expensive, it is often useful to estimate the number of irreducible \(\mathbf{q}\) points before running the phonon calculation (run a Gammapoint phonon calculation). Note however that this step is optional.
Step 3: nSCF calculation
We now run the nSCF calculations needed to generate the wavefunctions on the full \(\mathbf{k}\) point grid, which we’ll need both for generating Wannier functions with Wannier90 and for forming the coarsegrid eph matrix elements in Perturbo.
Here is the nscf.in input file (truncated):
&CONTROL
prefix = 'si'
calculation = 'nscf'
wf_collect = .true.
outdir='./tmp'
pseudo_dir='../pseudo'
/
&SYSTEM
ibrav = 2
celldm(1) = 10.264
nat = 2
nbnd=16
ntyp = 1
ecutwfc = 40.0
/
&ELECTRONS
conv_thr = 3.0d12
mixing_mode = 'plain'
mixing_beta = 0.7
diagonalization = 'david'
diago_full_acc = .true.
/
&IONS
/
&CELL
press_conv_thr = 0.01
/
ATOMIC_SPECIES
Si 28.085 Si_DOJO_PBEsol.upf
ATOMIC_POSITIONS crystal
Si 0.00000000 0.00000000 0.00000000
Si 0.25000000 0.75000000 0.25000000
K_POINTS crystal
512
0.00000000 0.00000000 0.00000000 1.953125e03
...
Make sure that the number of k points is commensurate with the number of \(\mathbf{q}\) points used for phonons, otherwise, qe2pert.x will stop. Remember to copy the QE ‘prefix’.save directory from the SCF calculation the current directory:
$ cp r ../scf/tmp ./
The k points lists could be produced by using kmesh.pl
in ‘utility’ folder of Wannier90:
$ kmesh.pl 8 8 8 >> nscf.in
Then run the nSCF calculation with QE (pw.x
executable).
Step 4: Wannier90 calculation
The directory contains two input files, one for wannier.x
and the other for pw2wannier90.x
.
Here is the pw2wan.in input file:
&inputpp
outdir='./tmp'
prefix = 'si'
seedname = 'si'
spin_component = 'none'
write_mmn = .true.
write_amn = .true.
write_unk = .false.
/
And here is the si.win input file (truncated):
begin projections
Si:sp3
end projections
guiding_centres=true
num_bands = 16
num_wann = 8
iprint = 2
dis_num_iter = 500
dis_win_min = 100.000
dis_win_max = 17.200
dis_froz_min = 100.000
dis_froz_max = 9.000
num_iter = 10000
mp_grid : 8 8 8
begin unit_cell_cart
bohr
5.1320 0.0000 5.1320
0.0000 5.1320 5.1320
5.1320 5.1320 0.0000
end unit_cell_cart
write_u_matrices = .true.
write_xyz = .true.
#restart = plot
BANDS_PLOT = TRUE
BANDS_PLOT_FORMAT = gnuplot
BANDS_NUM_POINTS = 100
BEGIN KPOINT_PATH
L 0.500 0.500 0.500 G 0.000 0.000 0.000
G 0.000 0.000 0.000 X 0.500 0.000 0.500
X 0.500 0.000 0.500 W 0.500 0.250 0.750
W 0.500 0.250 0.750 K 0.375 0.375 0.750
K 0.375 0.375 0.750 G 0.000 0.000 0.000
END KPOINT_PATH
begin atoms_frac
Si 0.00000 0.00000 0.00000
Si 0.25000 0.75000 0.25000
end atoms_frac
begin kpoints
0.0000000 0.0000000 0.0000000
...
end kpoints
In the input file si.win, we instruct Wannier90 to write two important quantities for qe2pert.x
, the \(U(\mathbf{k})\), \(U^{\text{dis}}(\mathbf{k})\) matrices and the position of the Wannier function centers, using: write_u_matrices=true
and write_xyz=true
.
We create tmp directory:
$ mkdir tmp
and change into it. We soft link to the QE ‘prefix’.save directory obtained in the nSCF calculation:
$ cd tmp
$ ln sf ../../nscf/tmp/si.save
We then run Wannier90. This step consists of three runs:

1) Run
wannier90.x pp si
to generate a list of the overlaps (the si.nnkp file, more details here). 
2) Run
pw2wannier90.x < pw2wan.in
. 
3) Run
wannier90.x si
(the si.win input file will be implicitly used).
The important output files for qe2pert.x
are si_u.mat, si_u_dis.mat, and si_centres.xyz. For disentangled bands, there would be no ‘prefix’_u_dis.mat. We encourage the user to become familiar with Wannier90 to run this step for different materials.
The user has to run Wannier90 3.0 or higher, since otherwise the \(U\) matrices will be printed out.
Step 5: Running qe2pert.x
We are ready to compute the eph matrix elements on the coarse \(\mathbf{k}\) point (determined by the nSCF step) and \(\mathbf{q}\) point (determined by the phonon step) Brillouin zone grids. First, copy or link the electronic and phonon calculation results to the current directory.
$ cd qe2pert
$ mkdir tmp
$ cd tmp
$ #link to the nscf .save directory
$ ln sf ../../pwphwann/nscf/tmp/si.save
$ cd ../
$ #link to the wannier information
$ ln sf ../pwphwann/wann/si_u.mat
$ ln sf ../pwphwann/wann/si_u_dis.mat
$ ln sf ../pwphwann/wann/si_centres.xyz
Here we show the input file (qe2pert.in) for the executable qe2pert.x
:
&qe2pert
prefix='si'
outdir='./tmp'
phdir='../pwphwann/phonon/save'
nk1=8, nk2=8, nk3=8
dft_band_min = 1
dft_band_max = 16
num_wann = 8
lwannier=.true.
load_ephmat = .false.
system_2d = .false.
/
The description of the input parameters:
 prefix: needs to be the same as the prefix used in the input files for QE.
 outdir: contains the save directroy obtained from the nSCF calculations. The calculated eph matrix elements will be stored in this directory.
 phdir: is the save directory inside which we collected all the phonon information.
 nk1, nk2, nk3: are the number of \(\mathbf{k}\) points along each direction used in the nSCF and Wannier90 calculations.
 dft_band_min and dft_band_max: determine the range of bands we are interested in, and should be the same as the values used in the Wannierization process. For example, if we used 40 bands in the nSCF calculation and we excluded bands 14 and 3140 in the Wannierization, then
dft_band_min=5
anddft_band_max=30
.  num_wann: the number of Wannier functions.
 lwannier: a logical flag. When it is
.true.
, the eph matrix elements are computed using the Bloch wave functions rotated with the Wannier unitary matrix; if.false.
, the eph matrix elements are computed using the Bloch wave functions, and the eph matrix elements are then rotated using the Wannier unitary matrix. By default, it is.true.
to reduce computational cost.  load_ephmat: a logical flag. If
.true.
, reuse eph matrix elements in Bloch function basis computed previously. This is useful if you want to test different Wannier bases. For example, you could first runqe2pert.x
withlwannier=.false.
, and then rerunqe2pert.x
withlwannier=.false.
andload_ephmat=.true.
with different Wannier unitary matrix.  system_2d: if the materials is twodimensional, so that in one direction only one \(\mathbf{k}\) point is used, set it to
.true.
; the default is.false.
.
Now we are ready to run the eph matrix elements:
export OMP_NUM_THREADS=4
$ mpirun np 2 qe2pert.x npools 2 i qe2pert.in > qe2pert.out
This task is usually timeconsuming on a single core, but it can be made much faster (minutes) on multiple cores.
The executables qe2pert.x
employ hybrid parallelization (MPI plus OpenMP), e.g. 2 MPI processes and each process span 4 OpenMP threads in this example.
To speed up the calculations, the users could increase the number of OpenMP threads and MPI processes.
Threads with OpenMP are particularly useful when the RAM (memory) of computing nodes is limited. The memory consumption reduces to minimum when using 1 MPI process per node and setting OMP_NUM_THREADS
to the number of cores per node.
Once the calcalculation has completed, we obtain the output file si_epr.h5, which is an HDF5 database with all the information needed to run perturbo.x
(which is described in the next section).
Structure of epr.h5 HDF5 file
The HDF5 file ‘prefix’_epr.h5 has the following data structure:
Variables  Data type  Meaning 

alat  real  lattice constant in atomic unit (Bohr) 
at  real; dimension (1:3,1:3)  lattice vectors in unit of \(a\), where \(a\) is the lattice constant 
bg  real; dimension (1:3,1:3)  reciprocal lattice vectors in unit of \(2\pi/a\), where \(a\) is the lattice constant 
nat  intger  number of atoms in the unit cell 
tau  real; dimension (1:3,1:nat)  atomic positions in a units, where \(a\) is the lattice constant 
volume  real  volume for the unit cell in unit of Bohr\(^3\) 
nsym  integer  number of symmetry operations for the simulation cell 
symop  integer; dimension(1:3,1:3,1:nsym)  symmetry operations 
kc_dim  integer; dimension (1:3)  number of coarse \(\mathbf{k}_c\)points in each direction for the electronic system 
spinor  integer  if spinor==1, the electron is spinor; if spinor==0, nonspinor 
polar_alpha  real  broadening of the Lorezian function for generating the weighting for polar electronphonon matrix element calculations 
epsil  real; dimesion (1:3,1:3)  the dielectric tensor 
qc_dim  integer; dimension (1:3)  number of coarse q\(_c\)points in each direction for the phononic system 
mass  real; dimension (1:nat)  atomic mass in atomic unit 
zstar  real; dimension (1:3,1:3,1:nat)  the effective charge for each atom 
system_2d  Integer  if system_2d==1, the system is a 2D material; if system_2d==0, not a 2D material 
num_wann  integer  number of Wannier functions used in the Wannierization 
wannier_center  real; dimesion (1:3,1:num_wan)  the center of Wannier functions in Cartesian unit 
wannier_center_cryst  real; dimesion (1:3, 1:num_wan)  the center of Wannier functions in crystal unit 
Variables  Datatype  Meaning 

hopping_rx  real; dimension(1:nr_el) where nr_el is the number of lattice points in the electronic system, and note that nr_el depends on which Wannier function is refered to  real part of the electronic Hamiltonian in Wannier basis; x is the index for the lattice points in the WignerSeitz cell 
hopping_ix  real; dimension(1:nr_el) where nr_el is the number of lattice points in the electronic system, and note that nr_el depends on which Wannier function is refered to  Imaginary part of the electronic Hamiltonian in Wannier basis; x is the index for the lattice points in the WignerSeitz cell 
Variables  Datatype  Meaning 

ep_hop_r_x_x_x  real; dimension (1:nat,1:nr_el,1:nr_ph) where nat is the number of atoms, nr_el is the number of lattice points in the electronic system, and nr_ph is the number of lattice points in the phononic system  real part of the electronphonon matrix elements in Wannier basis; 1st x: index for the atoms; 2nd x: index for R for electrons; 3rd x: index for R for phonons 
ep_hop_i_x_x_x  real; dimension (1:nat,1:nr_el,1:nr_ph) where nat is the number of atoms, nr_el is the number of lattice points in the electronic system, and nr_ph is the number of lattice points in the phononic system  imaginary part of the electronphonon matrix elements in Wannier basis; 1st x: index for the atoms; 2nd x: index for R for electrons; 3rd x: index for R for phonons 
Variables  Datatype  Meaning 

Ifcx  complex; dimesion (1:3,1:3,1:nr_ph) where nr_ph: number of lattice points for the phononic system  force constant 
The ‘prefix’_epr.h5 file structure can be schematically represented as follows: