The following calculation modes (`calc_mode`

parameter) are related to the e-ph scattering computation:

`'setup'`

: setup for transport calculations or carrier dynamics simulations.`'imsigma'`

: compute the e-ph self-energy for electronic crystal momenta read from a list.`'meanfp'`

: compute the e-ph mean free path, also output the corresponding band velocity and relaxation time.

## Setup electron k- and phonon q-grids

*calc_mode = ‘setup’*

**Computes:**Set up transport property calculations by providing \(\mathbf{k}\) points, \(\mathbf{k}\) point tetrahedra and (if needed) finding chemical potentials for given carrier concentrations. Also computes the density of states.

Requires to specify up to 14 variables in the input file (*pert.in*)

- prefix: same prefix as in
*‘prefix’_epwan.h5* - calc_mode: set to ‘setup’
- hole: By default,
`hole`

is set to`.false.`

. Set it to .true.**only**when computing hole mobility of a semiconductor. if*hole*is .true.,`perturbo.x`

computes hole concentration, instead of electron concentration. - boltz_kdim: number of \(\mathbf{k}\) points along each dimension of a \(\mathbf{k}\) point grid for the electrons momentum. This Gamma-centered Monkhorst-Pack \(\mathbf{k}\) point grid is employed to compute the mobility or conductivity.
- boltz_qdim: number of \(\mathbf{q}\) points along each dimension of a uniform grid for the phonon momentum; the default is that
`boltz_qdim(i)=boltz_kdim(i)`

. If users need the size as same as the \(\mathbf{k}\) grid, no need to specify these variables. Only phonons with mmentum on the \(\mathbf{q}\) grid are considered in the calculations of e-ph scattering. - boltz_emin, boltz_emax: energy window (in eV units) used to compute transport properties. The suggested values are from 6 k_BT below E_F (boltz_emin) to 6k_BT above E_F (boltz_emax), where E
_{F}is the Fermi energy, k_{B}the Boltzmann constant, and T is temperature in K units. - band_min, band_max: band window for transport property calculations
- ftemper: the filename of a file containing the temperature(s), chemical potential(s), and corresponding carrier concentration(s) for transport property calculations. Either chemical potentials or carrier concentrations is required dependending on the calculation setting.

Here is the input file (*pert.in*):

```
&perturbo
prefix = 'si'
calc_mode = 'setup'
boltz_kdim(1) = 80
boltz_kdim(2) = 80
boltz_kdim(3) = 80
boltz_emin = 6.4
boltz_emax = 6.9
band_min = 5
band_max = 6
ftemper = 'si.temper'
/
```

In the input file *pert.in*, we use a \(\mathbf{k}\) grid of 80 x 80 x 80 for electrons, which corresponds to boltz_kdim(i)=80, and use a \(\mathbf{q}\) grid for phonons of the same dimension as the \(\mathbf{k}\) grid. When a phonon \(\mathbf{q}\) grid different from the electron \(\mathbf{k}\) grid is desired, the user need to provide the \(\mathbf{q}\) grid variables *boltz_qdim(1)*, *boltz_qdim(2)*, and *boltz_qdim(3)* in the input file.

In this example, we want to compute the mobility of the electron carrier, so we choose an energy window that includes the conduction band minimum. Here the energy window is between 6.4 (*boltz_emin*) and 6.9 eV (*boltz_emax*), and the conduction band minimum is at 6.63 eV in this case. We include the two lowest conduction bands, with band indices 5 and 6 (*band_min* and *band_max*).

The ‘setup’ calculation finds all the relevant \(\mathbf{k}\) points (both irreducible and reducible \(\mathbf{k}\) points) and the tetrahedron needed for BZ integration for the given energy window and band window computes the DOS at the given energy window. It also computes carrier concentrations at given chemical potentials or determines the chemical potentials that corresponding to the given carrier concentrations, depending on the setting in the ftemper file.

In this case, the ftemper file `si.temper`

has the following format:

```
1 T
300.00 6.52 1.0E+18
```

The integer in the first line is the number of (temperature, chemical potential) settings at which we want to perform the transport calculations.
Each of the following lines contains three values, the temperature (K), Fermi level (eV), and carrier concentration (cm^{-3} in 3D materials or cm^{-2} in 2D materials).

The logical variable in the first line indicates whether to compute the carrier concentration for the input chemical potential (if `F`

) or determine the chemical potential corresponding to the input carrier concentration (if `T`

), thus only one of the chemical potential column and carrier concentration column in the *ftemper* file is meaningful.

The logical variable is only used in the `'setup'`

calculation. In all the other calc_mode options, `perturbo.x`

reads the chemical potential column and ignores the carrier concentration column (and the logical variable). If one wants to perform transport calculations at given carrier concentrations, then set the logical variable to `T`

in `'setup'`

calculations. `perturbo.x`

will find the corresponding chemical potentials and update the *ftemper* file accordingly (overwrite the chemical potential and carrier concentration columns and set the logical variable to `F`

).

**Note:**

`perturbo.x`

only search for chemical potentials within the given energy window, try extending the energy window if the updated *ftemper*file does not show reasonable carrier concentrations.

Run `perturbo.x`

with the following command (remember to link or copy *‘prefix’_epwan.h5* in the current directory):

```
$ mpirun -n 1 perturbo.x -npools 1 -i pert.in > pert.out
```

The calculation will take a few minutes or longer, depending the number of \(\mathbf{k}\) and \(\mathbf{q}\) points and the size of the energy window. We obtain 4 output files (*‘prefix’.doping*, *‘prefix’_tet.h5*, *‘prefix’_tet.kpt*, and *‘prefix’.dos*):

*‘prefix’.doping*contains chemical potentials and carrier concentrations for each tempearture of interest. The format is easy to understand so we do not show it here. Please take a look at the file by yourself.*‘prefix’_tet.h5*contains information on the \(\mathbf{k}\) points (both in the irreducible wedge and full grid) and the associated \(\mathbf{k}\) point tetrahedra in the energy window of interest. This file will be used to compute transport properties. Users familiar with HDF5 can read and manipulate this file with the standard HDF5 commands. The other users can just ignore the data stored in the file.*‘prefix’_tet.kpt*contains the coordinates (in crystal units) of the irreducible \(\mathbf{k}\) points in the energy window of interest. Note that the irreducible \(\mathbf{k}\) points coordinates is already included in*‘prefix’_tet.h5*, we output to this file in a format compatiable with that of*fklist*discussed in the calculation mode`'bands'`

(above) or`'imsigma'`

(below).*‘prefix’.dos*contains the density of states (number of states per eV per unit cell) as a function of energy (eV). The format is easy to understand so we do not show it here. The density of states sets the phase space for several electron scattering processes, so it is convenient to compute it and print it out.

In our example, since we used `'T'`

in the first line of ftemper, a new *ftemper* file is generated as output: that the *ftemper* file *‘si.temper’* has now become:

```
1 F
300.00 6.5504824219 0.9945847E+18
```

Note how `perturbo.x`

has computed the chemical potential (second entry in the second row) for the given temperature and carrier concentration (first and third entries of the second row). The logical variable in the first line is now `'F'`

, and *si.temper* can now be used as is in subsequent calculations.

The above explanation focuses on electrons. For holes carriers, please refer to *“example02-silicon-perturbo/perturbo/pert-setup-hole”*, link. In the input file for holes, remember to use `hole=.true.`

(default: `hole=.false.`

), and choose an appropriate energy window and the band indices for holes.

## Imaginary part of e-ph self-energy

*calc_mode = ‘imsigma’*

**Computes:**

The imaginary part of the lowest-order (so-called ‘Fan’) e-ph self-energy, \(\operatorname{Im}\Sigma\), for states in a range of bands and with crystal momenta \(\mathbf{k}\) read from a list (this list can be obtained from

`calc_mode='setup'`

or created manually). The scattering rates can also be obtained using \({2} \operatorname{Im}\Sigma /{\hbar}\). Variables in the input file (*pert.in*)

- prefix: the same prefix as the file
*‘prefix’_epwan.h5* - calc_mode: set to
`'imsigma'`

- band_min, band_max: bands used for transport property calculations
- ftemper: the filename of a file containing temperature, chemical potential, and carrier concentration values (see the format)
- fklist: the filename of a file containing the coordinates of a given electron \(\mathbf{k}\) point list (see the format)
- phfreq_cutoff: the cutoff energy for the phonons. Phonon with their energy smaller than the cutoff (in meV) is ignored; 0.5-2 meV is recommended.
- delta_smear: the broadening (in meV) used for the Gaussian function used to model the Dirac delta function
- fqlist: the filename of a file containing the coordinates of a given phonon \(\mathbf{q}\) point list will be used to compute the e-ph self-energy. For the format, see the section on the calculation mode
`'bands'`

. This is optional. If`fqlist`

is absent or`fqlist_=''`

, random \(\mathbf{q}\) points will be generated (see below). - sampling: sampling method for random \(\mathbf{q}\) points used in e-ph self-energy calculation. The default value is
`'uniform'`

, indicates sampling random \(\mathbf{q}\) points in the first BZ following uniform distribution. Another option is`'cauchy'`

, sampling random \(\mathbf{q}\) points following Cauchy distribution, which is useful for polar materials. Note that random \(\mathbf{q}\) points from other importance sampling methods or \(\mathbf{q}\) points on regular MP grid is also possible, one just needs to pre-generate the \(\mathbf{q}\) points list to a file, and pass the file to`perturbo.x`

via`fqlist`

. - cauchy_scale: the width of the Cauchy function; used only when sampling is ‘cauchy’.
- nsamples: number of random \(\mathbf{q}\) points sampled to compute the imaginary part of the e-ph self-energy for each \(\mathbf{k}\) point

Here is the input file (*pert.in*):

```
&perturbo
prefix = 'si'
calc_mode = 'imsigma'
fklist = 'si_tet.kpt'
ftemper = 'si.temper'
band_min = 5
band_max = 6
phfreq_cutoff = 1 ! meV
delta_smear = 10 ! meV
sampling = 'uniform'
nsamples = 1000000
/
```

In the current example, we compute the imaginary part of the e-ph self-energy of \(\mathbf{k}\) points in the *fklist* file (in this case, we use the irreducible Monkhorst-Pack \(\mathbf{k}\) point list in *si_tet.kpt* obtained from the calculation mode `'setup'`

). Note that if one is only interested in a high symmetry line, one can provide \(\mathbf{k}\) point path in the *fklist* file instead. The temperature, chemical potential for computing the e-ph self-energy are given in the *ftemper* file, *si.temper*, obtained from the perturbo `'setup'`

process (the carrier concentration column is ignored in `'imsigma'`

calculation). Note that `perturbo.x`

will do calculations, at once, for as many combinations of temperature and chemical potential as are specified in the lines below the first of *ftemper*.

Here we use a uniform random sampling (`sampling='uniform'`

) with 1 million random \(\mathbf{q}\) points (`nsample=1000000`

). The phonon frequency cutoff is 1 meV (`phfreq_cutoff=1`

), and the smearing for the Gaussian function is 10 meV (`delta_smear=10`

).

Before running `perturbo.x`

, remember to link or copy *‘prefix’_epwan.h5* in the current directory.

```
export OMP_NUM_THREADS=4
$ mpirun -n 8 perturbo.x -npools 8 -i pert.in > pert.out
```

This task is usually time-consuming on a single core, thus we run this calculation on multiple cores (32 cores in this case) using hybrid MPI plus openMP parallelization.

We obtain two output files:

*‘prefix’.imsigma*contains the computed imaginary part of the e-ph self-energy*‘prefix’.imsigma_mode*contains the computed imaginary part of the e-ph self-energy (where phonon modes are numbered for increasing energy values).

The following is the format of *‘prefix’.imsigma* (in this case, *si.imsigma*):

```
# Electron (Imaginary) Self-Energy in the Migdal Approx. #
# ( only for bands within [band_min, band_max] ) #
#----------------------------------------------------------#
# NO.k: 450 NO.bands: 2 NO.T: 1 NO.modes: 1
#
# Temperature(T)= 25.85203 meV; Chem.Pot.(mu)= 6.55048 eV
#===========================================================
# it ik ibnd E(ibnd)(eV) Im(Sigma)(meV)
#-----------------------------------------------------------
1 1 1 6.955370 1.2413716479777598E+01
......
......
#-----------------------------------------------------------
```

The variable *it* is a dummy variable for enumerating the temperature values, while, *ik* is the number of \(\mathbf{k}\) points in the fklist, *ibnd* the band number (in this case, band indices are 5 and 6). *Im(Sigma)* is the imaginary part of the e-ph self-energy (in meV units) for each state of interest.

Similarly, the format for *si.imsigma_mode* is

```
# Electron (Imaginary) Self-Energy in the Migdal Approx. #
# ( only for bands within [band_min, band_max] ) #
#----------------------------------------------------------#
# NO.k: 450 NO.bands: 2 NO.T: 1 NO.modes: 6
#
# Temperature(T)= 25.85203 meV; Chem.Pot.(mu)= 6.55048 eV
#===========================================================
# it ik ibnd E(ibnd)(eV) imode Im(Sigma)(meV)
#-----------------------------------------------------------
1 1 1 6.955370 1 1.415350400936959E+00
......
......
#-----------------------------------------------------------
```

Here we have an extra column with the phonon mode index (*imode*).

**Note:**One should always check the convergence of the e-ph self-energy with respect to the number of \(\mathbf{q}\) points and the smearing parameter (delta_smear). Check this paper for more detail.

Using the results in the *‘prefix’.imsigma* file, one can easily obtain, with a small script, the scattering rates for each state, which are equal to \({2}/{\hbar} \operatorname{Im}\Sigma\) (it’s convenient to use \(\hbar = 0.65821195\,\mathrm{eV}\,\mathrm{fs}\) to this end). Using additional tools provided in `perturbo.x`

, we can also compute the mean free path for each electronic state, as well as a range of phonon-limited transport properties.

One way of obtaining the relaxation times (and their inverse, the scattering rates) is to run the Python script `relaxation_time.py`

we provide to post-process the imsigma output (the desciption of the script is here). Another way is to obtain the relaxation times is to run a calculation of the mean free paths (see below), which conveniently outputs both the relaxation times and the mean free path for the desired electronic states.

Also note that an example calculation of the e-ph self-energy for holes, is provided in the example folder *“example02-silicon-perturbo/perturbo/pert-imsigma-hole”*, link, where we use different band indices (`band_min=2`

and `band_max=4`

), and the files, *fklist* and *ftemper*, are also different and obtained in a different perturbo `'setup'`

calculation.

## Electron mean free path

*calc_mode = ‘meanfp’*

**Computes:**The e-ph mean free paths for electronic states in a user-defined \(\mathbf{k}\) point list and range of bands.

**Note:**The mean free path calculation relies on the results of the calculation mode

`'imsigma'`

values obtained. Therefore, the user should first run the calculation mode `'imsigma'`

, and then compute the mean free pathsRequires the same files as `calc_mode='imsigma'`

but needs an additional file, *‘prefix’.imsigma*, obtained as an output in the `'imsigma'`

calculation.

Here is the input file (*pert.in*). It should be the same input as the one for the `'imsigma'`

calculation mode, except for the line specifying `calc_mode='meanfp'`

:

```
&perturbo
prefix = 'si'
calc_mode = 'meanfp'
fklist = 'si_tet.kpt'
ftemper = 'si.temper'
band_min = 5
band_max = 6
phfreq_cutoff = 1 ! meV
delta_smear = 10 ! meV
sampling = 'uniform'
nsamples = 1000000
/
```

Before running `perturbo.x`

, make sure you have the following files in the current directory (*“pert-meanfp-electron”*): *‘prefix’_epwan.h5*, *‘prefix’.imsigma* the *fklist* file (*si_tet.kpt* in this example), and the *ftemper* file (e.g., *si.temper* in this example). As explained above, one can reuse the input file of the calculation mode `'imsigma'`

by replacing the calculation mode with `calc_mode='meanfp'`

.

```
$ mpirun -n 1 perturbo.x -npools 1 -i pert.in > pert.out
```

This calculation usually takes only takes a few seconds. We obtain two output files:

*‘prefix’.mfp*contains the relaxation time and mean free path of each electronic state. Note that the MFP is the product of the state relaxation time and the absolute value of the band velocity.*‘prefix’.vel*contains the band velocity of each state

The format of *‘prefix’.mfp* is as follows:

```
#==========================================================#
# Electron Mean Free Path (tau_nk * |v_nk|, in nm) #
#==========================================================#
# NO.k: 2637 NO.bands: 2 NO.T: 1
#########
# Temperature(T)= 25.85203 meV; Chem.Pot.(mu)= 6.55048 eV
#-----------------------------------------------------------
# it ik ibnd E(ibnd)(eV) Relaxation time(in fs) MFP (in nm)
#-----------------------------------------------------------
1 1 1 6.955370 2.6511488462206518E+01 9.5929573542019302E+00
......
......
```

The variable *it* is the dummy variable for temperature; in this case, we only used one temperature (300 K). *ik* is the dummy variable for the given crystal momentum in the file fklist. *ibnd* is the dummy variable for bands; in this case, ibnd=1 corresponds to band index 5 and ibnd=2 is the band index 6. The 4^{th}, 5^{th}, and 6^{th} columns are energy (eV), relaxation time (fs), and mean free path (nm) of each state, respectively.

The format of *‘prefix’.vel* is shown below:

```
######################################################
# Band velocity #
######################################################
# ik ibnd E(ibnd)(eV) k.coord. (cart. alat) vel-dir |vel| (m/s)
1 1 6.955370 -0.01250 0.58750 -0.01250 -0.24926 -0.93581 -0.24926 3.6184152269976016E+05
......
......
```

The 1^{st} to 3^{rd} columns are the same as in *‘prefix’.mfp*. The 4^{th} to 6^{th} columns are the \(\mathbf{k}\) point coordinates in the crystal units. The 7^{th} to 9^{th} columns are the components of the unit vector specifying the direction of the velocity of each electronic states. The last column is the magnitude of the velocity (m/s) of each state.

For an example calculation of mean free paths for holes, please see the folder *“example02-silicon-perturbo/perturbo/pert-meanfp-hole”*, link.