Electron-phonon Scattering

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’

Directory: example02-silicon-perturbo/perturbo/pert-setup-electron/ link

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’_epr.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 momentum 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 6k_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.
find_efermi: whether or not to find the chemical potential using the carrier concentrations inputted via the ftemper file.

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

 find_efermi = .true. 
 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. We also set find_efermi to .true., indicating that we want the calculation to find the chemical potential using the inputted carrier concentration (the carrier concentration will be inputted using the ftemper file described below).

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 eV (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 Brillouin zone integration for the given energy window and band window computes the DOS at the given energy window.

In this case, the ftemper file si.temper has the following format:

1
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).

Note that depending on the value of find_efermi from pert.in, the calculation will either compute the carrier concentration for the input chemical potential (if .false.) or determine the chemical potential corresponding to the input carrier concentration (if .true.), thus only one of the chemical potential column and carrier concentration column in the ftemper file is meaningful. Here, we have set find_efermi to ‘.true.’

The parameter find_efermi 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. If one wants to perform transport calculations at given carrier concentrations, then set find_efermi to .true. 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).

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’_epr.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 the YAML file, ‘prefix’_setup.yml and the HDF5 file ‘prefix’_tet.h5.

The ‘prefix’_tet.h5 file 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 outputted YAML file, here called si_setup.yml, contains the inputs and outputs of the 'setup' calculation. The fields grid, DOS, and carrier density contain the information unique to the 'setup' calculation. The remaining fields are common to all Perturbo calculations, and their format is described here.

The grid, DOS, and carrier density fields are given below.

grid:
   number of  tetrahedra selected:      76296
   number of reducible k-points:      15975
   number of irreducible k-points:        450


DOS:

   energy units: eV

   DOS units: num. of states/eV/u.c.

   energy:
      -   0.663261E+01
      ......
      -   0.690061E+01

   DOS:
      -     0.0000000000E+00
      ......
      -     0.2673882268E+00

carrier density:

   temperature units: K

   chemical potential units: eV

   concentration units: cm-3

   # Cofiguration means one line in the temper file:
   # temperature efermi concentration
   number of configurations:   1

   configuration index:

      1:
         temperature: 300.00
         chemical potential:   6.5504824219
         concentration:  0.9945835E+18

The following information is given:

The grid field gives information on the number of tetrahedra, reducible \(\mathbf{k}\) points, and irreducible \(\mathbf{k}\) points in the energy window of interest. The specific \(\mathbf{k}\) points and tetrahedra are given in the ‘prefix’_tet.h5 file.
Under the DOS field, the energy field gives the energies (in eV units) for which the density of states will be computed, and the inner DOS field gives the corresponding density of states (number of states per eV per unit cell) as a function of those energies. The density of states sets the phase space for several electron scattering processes, so it is convenient to compute it and print it out.
Under the carrier density field, information is given on the temperatures (K units), chemical potentials (eV units) and carrier concentrations (cm^-3 units) for each of the configurations defined from the ftemper file.

In our example, since we set find_efermi to .true. in the pert.in file, a new ftemper file is generated as output: that the ftemper file ‘si.temper’ has now become:

1
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 file si.temper can now be used as is in subsequent calculations.

Click to show legacy-format ASCII text outputs of the setup calculation

Alongside the YAML and HDF5 files, we generate three additional output files ('prefix.doping,' 'prefix_tet.kpt,' and 'prefix.dos'), each containing a duplicate set of legacy-format information.

'prefix'.doping contains chemical potentials and carrier concentrations for each temperature of interest. The format is easy to understand so we do not show it here.
'prefix'_tet.kpt contains the coordinates (in crystal units) of the irreducible k points in the energy window of interest.
'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 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’

Directory: example02-silicon-perturbo/perturbo/pert-imsigma-electron/ link

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’_epr.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’_epr.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.

The outputted YAML file, here called si_imsigma.yml, contains the inputs and outputs of the 'imsigma' calculation. The field labeled imsigma contains the information unique to the 'imsigma' calculation. The remaining sections are common to all Perturbo calculations, and their format is described here.

The first part of the imsigma field is given below.

imsigma:

   temperature units: meV

   energy units: eV

   Im(Sigma) units: meV

   chemical potential units: eV

   k-point coordinate units: crystal

   number of k-points:     450

   # Band count starts from band_min
   number of bands:   2

   number of phonon modes:   6

   # Cofiguration means one line in the temper file:
   # temperature efermi concentration
   number of configurations:   1

   k-point coordinates:
      - [   0.00000,     0.28750,     0.30000,  ]
      ......
      - [   0.06250,     0.45000,    -0.48750,  ]

   energy:

      band index:

         1:
            -    6.9553696604
            ......
            -    7.0300563001
         2:
            ......

The first several fields of the imsigma section contain information on the units and inputs to the calculation. The three-dimensional \(\mathbf{k}\) points are given in crystal coordinates under the k-point coordinates field (each row is a point). Under the energy field and then band index field, there are two band indices, one for each of the two bands used in the calculation (i.e. bands between band_min and band_max). Under each band index, the band energies at each of the \(\mathbf{k}\) points are given in eV units. Note that band index 1 corresponds to 5^th band of silicon, and band index 2 corresponds to the 6^th band of silicon since only bands 5 and 6 were used in the calculation.

The main results of the 'imsigma' calculation are given in the configuration index field under the imsigma field, shown below.

   configuration index:

      1:

         temperature:  25.85203

         chemical potential:   6.5504824219

         band index:

            1:

               Im(Sigma):

                  total:
                     -  1.2224418754808822E+01
                     ......
                     -  1.4777769003852226E+01

                  phonon mode:

                     1:
                        -  1.3602930790051577E+00
                        ......
                        -  1.8789546796477412E+00
                     ......
                     6:
                        -  5.4880944318394089E+00
                        ......
                        -  3.2125090317485521E+00

            2:
               ......

In the configuration index field, there is only one field entitled ‘1’ because we only used one configuration. The following information is given for each of the configurations defined from the ftemper file.

The temperature (K units) and chemical potential (eV units).
Under the band index field and then the Im(Sigma) field, the imaginary part of the e-ph self-energy is listed for each of the bands between band_min and band_max (here there are 2 bands). These values are computed in total (under the total field) and for each of the 6 phonon modes (under the phonon mode field). Note that phonon modes are numbered for increasing energy values. Again, band index 1 corresponds to 5^th band of silicon, and band index 2 corresponds to the 6^th band of silicon.

We may next use Perturbopy to export the data from the YAML file to Python for postprocessing. For more details, see the Perturbopy tutorial.

Click to show legacy-format ASCII text outputs of the imsigma calculation

We obtain two additional ASCII files containing a copy of the ``'imsigma'`` calculation results.

'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 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.

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’

Directory: example02-silicon-perturbo/perturbo/pert-meanfp-electron/ link

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 paths

Requires 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’_epr.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

The outputted YAML file, here called si_meanfp.yml, contains the inputs and outputs of the 'meanfp' calculation. The field labeled “meanfp” contains the information unique to the 'meanfp' calculation. The remaining sections are common to all Perturbo calculations, and their format is described here.

The first part of the meanfp field is given below.

meanfp:

   temperature units: meV

   chemical potential units: eV

   energy units: eV

   MFP units: nm

   band velocity units: m/s

   relaxation time units: fs

   k-point coordinate units: crystal

   number of k-points:     450

   # Band count starts from band_min
   number of bands:   2

   # Cofiguration means one line in the temper file:
   # temperature efermi concentration
   number of configurations:   1

   k-point coordinates:
      - [  -0.01250,     0.58750,    -0.01250,  ]
      ......
      - [   0.87500,    -0.10000,     1.00000,  ]

   energy:

      band index:

         1:
            -    6.9553696604
            -    7.0300563001

         2:
            ......

The first part of the meanfp field contains information on the units and inputs to the calculation. The three-dimensional \(\mathbf{k}\) points are given in crystal coordinates (each row is a point). Under the energy and then band index field, the band energies at each of these \(\mathbf{k}\) points are given in eV units for bands between band_min and band_max. Here, there are 2 band indices, reflecting the 2 bands used in the calculation. Note that band index 1 corresponds to 5^th band of silicon, and band index 2 corresponds to the 6^th band of silicon since only those two bands were used in the calculation.

The next fields under the meanfp field are given below.

   band velocity:

      band index:

         1:
            -  3.6184152269975992E+05
            ......
            -  7.4153127777954191E+05
 
         2:
            ......

   band velocity direction:

      band index:

         1:
             - [  -0.24926,    -0.93581,    -0.24926,  ]
             ......
             - [  -0.99800,     0.06323,     0.00000,  ]

         2:
             ......

Under the band velocity field and then the band index field, the magnitude of the band velocities at each electronic state are given in m/s units for bands between band_min and band_max. The band velocity direction field is structured similarly and gives the components of the unit vector specifying the direction of the velocity of each electronic states.

The remaining results are given under the configuration index field below. Note that in this case, we only used one configuration.

   configuration index:

      1:

         temperature:  25.85203

         chemical potential:   6.5504824219

         band index:

            1:

               MFP:
                     -  9.7415063395807646E+00
                     ......
                     -  1.6588387141733620E+01

               relaxation time:
                     -  2.6922024500941074E+01
                     ......
                     -  2.2370448339557921E+01

            2:
               ......

Under the configuration index field, data for each of the configurations from the ftemper file are given:

The temperature of the configuration is given in K units, and chemical potential is given in eV units.
Under the band index field, and then the MFP field, the mean free paths for the states at each of the \(\mathbf{k}\) points are given in nm units for each band between band_min and band_max. Note that the MFP is the product of the state relaxation time and the absolute value of the band velocity.
Under the band index field, and then the relaxation time field, the relaxation times for the states at each of the \(\mathbf{k}\) points are given in fs units for each band between band_min and band_max.

We may next use Perturbopy to export the data from the YAML file to Python for postprocessing. For more details, see the Perturbopy tutorial.

Click to show legacy-format ASCII text outputs of the meanfp calculation

We obtain two output files containing a copy of the results:

'prefix'.mfp contains the relaxation time and mean free path of each electronic state.
'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 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.

Setup electron k- and phonon q-gridscalc_mode = ‘setup’

Imaginary part of e-ph self-energycalc_mode = ‘imsigma’

Electron mean free pathcalc_mode = ‘meanfp’

Setup electron k- and phonon q-grids
calc_mode = ‘setup’

Imaginary part of e-ph self-energy
calc_mode = ‘imsigma’

Electron mean free path
calc_mode = ‘meanfp’