The interface between Quantum Espresso and perturbo-ed: qe2pert_ed

Preparation for qe2pert_ed calculations

The essential ingredients for the interface qe2pert_ed are similiar to those for the interface qe2pert, except that the supercell calculations are needed in advance. Here we would recommend users to create additional directory called supercell in the pw-ph-wan direcotry. In the following we would use two types of defects in silicon–one is neutral silicon and the other one is charged phosphorus–as examples to demonstrate how to use the interface.

In the supercell directory, we create three directories: 1) pristine; 2) neutral-vacancy; 3) charged-P-doped. To reduce the computational cost, we use a supercell size of \(4\times 4\times 4\); note that the number means the supercell is constructed by repeating the primitive cell along the direction of each lattice vector. In the pristine directory, we create a directory called 040404 to indicate the size of the supercell. In the neutral-vacancy directory, we create two directories: 1) 040404-fixed; 2) 040404-relaxed, which mean the atoms around the defect are fixed or relaxed, i.e., neutral vacancy in this case. We do a similar thing for the charged-P-doped directory. The reason why we need the “fixed” case is to compute the potential alignment while the “relaxed” one is used to construct the electron-defect perturbation potential. Those supercell calculations can be done using the Gamma point to save the computational cost.

Users can find the above supercell example input files here: example08-silicon-defects/pw-ph-wan/supercell. To demonstrate the example, we use the same input files for the fixed and the relaxed case.

Calculations modes for qe2pert_ed

Once all the scf, nscf, phonon, wannierization, and supercell calculations for silicon are ready, we can proceed to run the interface between QE and electron-defect perturbo, i.e., qe2pert_ed.x, which generates an HDF5 file, called prefix_edwan.h5, used to run transport calculations in the main perturbo codes.

The interface has the following calculation modes (calc_mode):

1. calc_mode = ‘vqbase’, which is used to generate the interpolation table for the neutral-part of the defect perturbation potential (an expensive step)

2. calc_mode = ‘edmel’, which is used to compute the e-d matrix elements in the Bloch basis on a coarse k grid (an expensive step)

3. calc_mode = ‘edwan’, which is used to convert the e-d matrix elements from the Bloch representation to the Wannier representation

4. calc_mode = ‘inv_scr’, which is used only for charged defects to compute the inverse screening length for an isolated point charged defect; the current version only support cubic systems.

5. calc_mode = ‘edbloch_line’, which is used to directly compute the matrix elements using the wave functions along the given high symmetry lines that specified in the nscf calculations.

In the following, we look into the neutral vacancy and the charged phophorus in silicon, separately, for holes, which are computationally cheaper to demonstrate. First, we create a directory called qe2perted, at the same directory level as the pw-ph-wann directory, to put all the qe2perted caculation input and output files. For each defect type, we would recommend to create a corresponding directory inside the qe2perted directory; in this case, neutral-vacancy-040404 and charged-P-doped-040404. Inside each defect directory, please link the wannier unitary matrices and wannier centers via soft-link as follows.

Inside the qe2perted/neutral-vacancy-040404 directory (similarly for the charged defect):

>> ln -sf ../../pw-ph-wann/wann/si_u.mat ./
>> ln -sf ../../pw-ph-wann/wann/si_centres.xyz ./

Then we create another directory called tmp in the defect directory. In the tmp directory, we link the prefix.save in the nscf calculation as well as the supercell calculations:

>> ln -sf ../../../pw-ph-wann/nscf/tmp/si.save ./
>> ln -sf ../../../pw-ph-wann/supercell/pristine/040404/si-u040404.save ./
>> ln -sf ../../../pw-ph-wann/supercell/neutral-vacancy/040404-fixed/tmp/si-u040404-neutral-vacancy-fixed.save ./
>> ln -sf ../../../pw-ph-wann/supercell/neutral-vacancy/040404-relaxed/tmp/si-u040404-neutral-vacancy-relaxed.save ./

Once we are done with the proper links, we can proceed to run qe2perted calculations. The input file for qe2perted calculations has a namelist called qe2pert_ed, as shown in the following input file examples.

Neutral vacancy in silicon

The first step is to run the calculation mode: “vqbase”

An input file example (“pert-ed-vqbase.in”) for the calculation mode “vqbase” is as follows. The meaning of the variables are explained in the comment lines. The input file can be found here: example08-silicon-defects/qe2perted/neutral-vacancy-040404

&qe2pert_ed
  calc_mode = 'vqbase'

  ! output directory where the "QE .save" directories for
  ! the primitive and the supercells are linked
  outdir = './tmp'

  ! the prefix for the primitive cell (pmcell)
  pmcell = 'si'

  ! the pristine supercell 
  pscell = 'si-u040404'
  ! supercell containing the fix defect 
  ! i.e., without the atomic relaxation
  dscell_fix = 'si-u040404-neutral-vacancy-fixed'
  ! supercell containing the relaxed defect
  dscell_rel = 'si-u040404-neutral-vacancy-relaxed'

  ! the defect position (unit: angstrom)
  dcenter(1) = -10.86199805
  dcenter(2) =  10.86199805
  dcenter(3) =  10.86199805

  ! potential alignment method
  ! 1: core-average; 2: voronoi-cell
  align_num = 1

  ! the range of q for V(q), the fourier coefficients for the 
  ! defect perturbation potential; here q is in unit of Bohr^{-1} 
  qrange(1) = 6.0
  qrange(2) = 6.0
  qrange(3) = 6.0
  ! the q step 
  qstep = 0.1
/

The defect position must be provided in unit of angstrom. The variable for the potential alignment calculation for the supercells, align_num, can be computed using the core-average (align_num =1, the default) or the voronoi-cell (align_num = 2) method. The variable qrange is the range of the q vectors (along the x, y, and z direction) for V(q), which is the Fourier coefficients of the defect perturbation potential, while the variable qstep is the step size along each direction within the given range. The suggested values for the qrange and qstep are $2\pi$ and $0.1$, respectively. To check whether the values for the two variables, please refer to the section for the calculation mode “vqcheck”.

After the calculation, the code gives an HDF5 file, “si_edwan.h5”, which contains only the information of the defect perturbation potential. We still need to run the other calculation modes to get additional information before running the perturbo transport calculations.

The second step is to run the calculation mode: “edmel”

An input file example (“pert-ed-edmel.in”) for the calculation mode “edmel” is as follows. The meaning of the variables are explained in the comment lines. The input file can be found here: example08-silicon-defects/qe2perted/neutral-vacancy-040404

&qe2pert_ed
  calc_mode = 'edmel'
  outdir = './tmp'

  pmcell = 'si'

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-neutral-vacancy-fixed' 
  dscell_rel = 'si-u040404-neutral-vacancy-relaxed'

  dcenter(1) = -10.86199805
  dcenter(2) =  10.86199805
  dcenter(3) =  10.86199805

  ! provide the information for e-d matrix elements calculations 
  ! for exmaple, how many k points on the coarse grid, whether to 
  ! use the wannier gauge in advance, the band index, and number of 
  ! wannier functions 

  ! k point along each direction used in the wannierization 
  nk1 = 8, nk2 = 8, nk3 = 8

  ! band index used in the DFT 
  dft_band_min = 1,
  dft_band_max = 4,

  ! number of wannier functions 
  num_wann = 4

  ! wannier gauge
  lwannier = .true.

  ! method to compute the plane-wave matrix
  ! plm_num: 1 => compute in Fourier space directly
  ! plm_num: 2 => compute in real space first and then convert back to k space (slower)
  plm_num = 1

  ! the range of G vector used in the e-d matrix element calculations
  ! in principles, it should be infinity, but in practice, it only needs 
  ! a short range; users need to find the optimal condition (unit: Borh^-1) 
  ! the way to find the optimal condition is to use
  ! another calculation mode "edbloch_line"
  gmcut = 6.0
/

The additional variables, compared to the previous input file for the calculation mode “vqbase”, include the information used in the wannierization process, for example, the number of k points in the coarse k grid (i.e., nk1, nk2, and nk3) and the number of wannier functions (i.e., num_wann). The variable lwannier (the default is .true.) means the Wannier gauge is applied to the e-d matrix element calculations before the Wannier transformation, and it can reduce the computational cost by using less DFT electronic states.

Computing the e-d matrix elements on the k coarse grid involves the calculations of the plane-wave matrix elements, which can be computed using two approaches: one (plm_num = 1, the default and the faster) uses the Fourier coefficients of the wave functions on the coarse grid, while the other (plm_num = 2) uses the wave functions in real space, which is then converted into the k space. The last variable gmcut is used to constrain the summation of the G vectors used in the e-d matrix element calculations; in principles, it should be the same as the set of the G vectors used in the scf calculations, but in practice, the results can be converged within a small set of the G vectors.

The matrix elements computed on the coarse grid are stored in the file called prefix_eld.h5 in the outdir directory, which is specify in the input file; in this case, “tmp/si_eld.h5”, which will be read when we run the next calculation mode “edwan”.

The last step is to run the calculation mode: “edwan”

An input example (“pert-ed-edwan.in”) for the calculation mode “edwan”: The input file can be found here: example08-silicon-defects/qe2perted/neutral-vacancy-040404

&qe2pert_ed
  calc_mode = 'edwan'
  outdir = './tmp'

  pmcell = 'si'

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-neutral-vacancy-fixed'
  dscell_rel = 'si-u040404-neutral-vacancy-relaxed'

  dcenter(1) = -10.86199805
  dcenter(2) =  10.86199805
  dcenter(3) =  10.86199805

  ! specify the size of the supercell
  nsu(1) = 4, nsu(2) = 4, nsu(3) = 4

  nk1 = 8, nk2 = 8, nk3 = 8

  dft_band_min = 1,
  dft_band_max = 4,
  num_wann = 4
  lwannier = .true.

  plm_num = 1
  gmcut = 6.0
/

Compared to the input file for the calculation mode “edmel”, an additional information needs to be provided is the size of the supercell, specified by the variable nsu. This calculation mode runs faster than the previous two ones.

After the calculation, the e-d matrix elements in the wannier basis are stored in the HDF5 file prefix_edwan.h5, in this case, si_edwan.h5, which is used in the perturbo calculations.

Additional checking: to check the interpolation table by running the calculation mode: ‘vqcheck’, for both neutral and charged defects

The purpose of this calculation mode is used to check where the variables qrange and qstep for constructing the interpolation table of the Fourier coefficients, V(q), are good enough to reproduce those directly computed ones. Here we provide an input file example (“pert-ed-vqcheck.in”) as follows. The input file can be found in example08-silicon-defects/qe2perted/neutral-vacancy-040404 and example08-silicon-defects/qe2perted/neutral-vacancy-040404.

&qe2pert_ed
  calc_mode = 'vqcheck'
  outdir = './tmp'
  pmcell = 'si'

  fvqlist = 'qpath.dat'
/

Here the varible fvqlist points to the file ‘qpath.dat’ that contains the q points used for computing V(q). The format of the file is similar to that used in the perturbo “ephmat” calculations (please see the link: XXX). Note that the q points shoule be within the qrange specified in the “vqbase” calculation mode.The output file prefix.vqline (in this case, si.vqline) contains the interpolated and the directly computed V(q).

Additional checking: to check e-d matrix elements by running the calculation mode: ‘edbloch_line’, for both neutral and charged defects

The purpose is to compute the e-d matrix elements in the Bloch basis, which is used as reference for the interpolated e-d matrix elements as well as to find the optimal value for the variable gmcut. This calculation must be performed after the mode “vqbase”, because it requires the perturbation potential.

Here we need to run a nscf calculation with a QE input file, which has an initial k point followed by a few final k points that we are interested in. Therefore, we create a directory called nscf-path in the pw-ph-wann directory, and run the nscf calculation there. Here is the input file example: example08-silicon-defects/pw-ph-wan/nscf-path. Once we have the nscf results, inside the defect directory in the qe2perted directory, we create a directory called tmp-path, and provide the soft-link to the .save directory of the nscf-path directory and the supercell .save directories. For example, inside the tmp-path directory, use the following command:

>> ln -sf ../../../pw-ph-wann/nscf-path/tmp/si.save ./
>> ln -sf ../../../pw-ph-wann/supercell/pristine/040404/si-u040404.save ./
>> ln -sf ../../../pw-ph-wann/supercell/neutral-vacancy/040404-fixed/tmp/si-u040404-neutral-vacancy-fixed.save ./
>> ln -sf ../../../pw-ph-wann/supercell/neutral-vacancy/040404-relaxed/tmp/si-u040404-neutral-vacancy-relaxed.save ./

Then we are ready to compute the corresponding matrix elements. The input file for the calculation mode can be found here: example08-silicon-defects/qe2perted/neutral-vacancy-040404

&qe2pert_ed
  calc_mode = 'edbloch_line'
  outdir = './tmp-path'
  pmcell = 'si'

  nk1 = 8, nk2 = 8, nk3 = 8
  dft_band_min = 1,
  dft_band_max = 4,
  num_wann = 4
  lwannier = .false.

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-neutral-vacancy-fixed'
  dscell_rel = 'si-u040404-neutral-vacancy-relaxed'

  dcenter(1) = -10.86199805
  dcenter(2) =  10.86199805
  dcenter(3) =  10.86199805

  plm_num = 1
  gmcut = 6.0

  band_min = 1
  band_max = 4
/

After the calculation, it gives a file prefix.edbline (in this case, si.edbline), which contains the absolute value of the matrix elements that sum over the bands specified by the variables band_min and band_max.

One can use this mode: 1) to check what value of the variable gmcut is optimal to converge the matrix elements, and 2) to compare the interpolated matrix elements obtained later in the perturbo calculations.

charged phophorous in silicon

The first step is to run the calculation mode: “inv_scr”

Here we compute the inverse of the screening lengths from the density of states of the primitive cell, given the doping concentration(s) and the temperature(s) in the file called si.temper, which is specified using the variable ftemper in the input file; the format for the file is the same as that used in the perturbo setup calculation. The example ftemper file can be found: example08-silicon-defects/qe2perted/charged-P-doped-040404

An input file example (“pert-ed-inv-scr.in”) for the calculation mode “inv_scr” is as follows. The input file can be found here: example08-silicon-defects/qe2perted/charged-P-doped-040404

&qe2pert_ed
  calc_mode = 'inv_scr'
  outdir = './tmp'
  pmcell = 'si'

  nk1 = 8, nk2 = 8, nk3 = 8

  dft_band_min = 1,
  dft_band_max = 4,
  num_wann = 4
  lwannier = .true.

  boltz_emin = 6.0
  boltz_emax = 6.6

  boltz_kdim(1) = 160
  boltz_kdim(2) = 160
  boltz_kdim(3) = 160

  hole = .true.

  ftemper = 'si.temper'

  ! static dielectric constant (only for cubic systems)
  epsc = 11.3

  band_min = 2
  band_max = 4

  ! if consider incomplete ionization
  ! turn on ldop
  !ldop = .true.
  !!VBM = 6.2286
  !dlev = 6.2656
/

Here we need to provide the static dielectric constant for the variable epsc, which means the epsilon constant \(\varepsilon(0)\). To compute the density of states, the energy window is specified by the two variables boltz_emin and boltz_emax, and the number of k points in the fine grid is specified by the varibale boltz_kdim.

The calculation gives an output file called prefix.inv_scr_length (in this case, si.inv_scr_length), where the inverse of the screening length is stored and used later for the e-d matrix element calculations.

By default, we assume that all the dopants are completely ionized. If we consider the incompletely ionized case, the two additional variables, ldop and dlev, needs to be specified. When ldop is .true., it means the code considers the incompletely ionized case. The variable dlev means the energy of the defect level; here it is computed using the conduction band minimum minus the ionization energy of the phophorus dopants (around 45 meV) from literature.

The following steps are the same as those for the neutral defects. Also, the input file variables are very similar, except those used for the information of the charged defect.

The second step is to run the calculation mode: “vqbase”

An input example (“pert-ed-vqbase.in”) for the calculation mode “vqbase” is as follows. An additional variables for the charged defect information, compared to the neutral defect case, needs to provided. The input file can be found here: example08-silicon-defects/qe2perted/charged-P-doped-040404

&qe2pert_ed
  calc_mode = 'vqbase'
  outdir = './tmp'
  pmcell = 'si'

  ! make sure your system has the cubic symmetries.
  ! the default is .false. => the code terminates
  lcube =.true.
  ! charged defect information
  ! the defect is a charged defect 
  lcharged = .true.
  ! the charge value of the defect, which must the same
  ! as the one used in the supercell calculation 
  zcharge = 1.0
  ! static dielectric constant (only for cubic systems)
  epsc = 11.3

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-P-charged-fix'
  dscell_rel = 'si-u040404-P-charged-fix'

  dcenter(1) = -6.78875
  dcenter(2) =  6.78875
  dcenter(3) =  6.78875

  align_num = 1

  qrange(1) = 6.0
  qrange(2) = 6.0
  qrange(3) = 6.0
  qstep = 0.1
/

The variables for the charged defects include: 1) lcharged tells qe2pert_ed.x that we are dealing with a charged defect; 2) zcharge should be the same as the variablse tot_charge used in the QE input file for the supercell calculation; 3) epsc means the static dielectric constant (currently, only for cubic systems). Note that the variable lcube must be turned on to make sure that users know they are running a cubic system for charged defects.

The calculation stores the interpolation table for the neutral defect perturbation potential in the HDF5 file, in this case, si_edwan.h5.

The third step is to run the calculation mode: “edmel”

An input example (“pert-ed-edmel.in”) for the calculation mode “edmel” is as follows. The input variables are similar to the neutral defect case, except those variables for the charged defect, which have been described above. The output files are similar to the neutral defect case. The input file can be found here: example08-silicon-defects/qe2perted/charged-P-doped-040404

&qe2pert_ed
  calc_mode = 'edmel'
  outdir = './tmp'
  pmcell = 'si'

  lcube =.true.
  lcharged = .true.
  zcharge = 1.0
  epsc = 11.3

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-charged-P-fixed'
  dscell_rel = 'si-u040404-charged-P-relaxed'

  dcenter(1) = -6.78875
  dcenter(2) =  6.78875
  dcenter(3) =  6.78875

  nk1 = 8, nk2 = 8, nk3 = 8

  dft_band_min = 1,
  dft_band_max = 4,
  num_wann = 4
  lwannier = .true.

  plm_num = 1
  gmcut = 6.0
/

The forth step is to run the calculation mode: “edwan”

Here the input file is basically the same as that for neutral defects, except the variables for the charged defect information, described above. An input example (“pert-ed-edwan.in”) for the calculation mode “edwan” is as follows. Once this step is done, we are ready to compute the transport properties due to the e-d interaction from charged defects in perturbo calculations. The input file can be found here: example08-silicon-defects/qe2perted/charged-P-doped-040404

&qe2pert_ed
  calc_mode = 'edwan'
  outdir = './tmp'
  pmcell = 'si'

  lcube =.true.
  lcharged = .true.
  zcharge = 1.0
  epsc = 11.3

  pscell = 'si-u040404'
  dscell_fix = 'si-u040404-charged-P-fixed'
  dscell_rel = 'si-u040404-charged-P-relaxed'

  nsu(1) = 4, nsu(2) = 4, nsu(3) = 4

  dcenter(1) = -6.78875
  dcenter(2) =  6.78875
  dcenter(3) =  6.78875

  nk1 = 8, nk2 = 8, nk3 = 8

  dft_band_min = 1,
  dft_band_max = 4,
  num_wann = 4
  lwannier = .true.

  plm_num = 1
  gmcut = 6.0
/

Peturbo for electron-defect interaction

Once we have the HDF5 prefix_edwan.h5 obtained from the interface qe2pert_ed.x, we can proceed to perform the perturbo transport calculations. Users that are familiar with the perturbo caluclations for e-ph calculations do not have to learn any new variables, except two new variables leph and led, which are the logic flags to specify whether to run the e-ph, the e-d or both calculations. The default value for leph is .true., which requires the e-ph HDF5 file, prefix_epwan.h5, while the default one for led is .false., which requires the e-d HD5F file, prefix_edwan.h5. Additional calculation mode is “edmat”, which is analogous to “ephmat” in the e-ph matrix element calculations.

Remember that to run the perturbo calculations, make sure using the soft link to link to the corresponding HDF5 file.

Here the perturbo input files for neutral and charged defects are basically the same, except the calculation modes that require the ftemper file. For charged defects, the ftemper file must be the same as that used in the qe2pert_ed calculations. Therefore, we will look into each perturbo calculation mode for the two types of defects together. First, we create two directories called perturbo-netural-vacancy-040404 and perturbo-charged-P-doped-040404 for the neutral defect and the charged defect, respectively. In each directory, we suggest users use the calculation mode with a prefix “pert-“ as the names for the subdirectories. For example, for the calculation mode “ calc_mode = ‘edmat’ “, we create a directory called “pert-edmat”.

The example input files for the following calculation modes can be found here: 1) example08-silicon-defects/perturbo-netural-vacancy-040404 and 2) example08-silicon-defects/perturbo-charged-P-doped-040404

calculation mode: ‘edmat’

The purpose is to compute the interpolated e-d matrix elements for the given k and k’ points, which are specified in two files with the input variable names fklist and fqlist, the format of which are the same as that used in the e-ph perturbo calculation mode “ephmat”.

In our example here, the two files are si_klist.kpt and si_kplist.kpt. Here we use one k point for the initial state and a list of several k’ points along the high symmetry lines for the finial states.

An input file example (“pert.in”) for the neutral defect (in this case, netural vacancy in silicon) looks as follows:

! for netural defect 
&perturbo
   calc_mode = 'edmat'
   prefix = 'si'

   leph = .false.
   led = .true.

   band_min = 1
   band_max = 4

   fklist = 'si_klist.kpt'
   fqlist = 'si_kplist.kpt'
/

An input file example (“pert.in”) for the charged defects (in this case, phosphorous dopant in silicon) looks as follows:

! for charged defect
&perturbo
   calc_mode = 'edmat'
   prefix = 'si'

   leph = .false.
   led = .true.

   band_min = 1
   band_max = 4

   fklist = 'si_klist.kpt'
   fqlist = 'si_kplist.kpt'

   ! needed for charged defects, and the file must be the same as 
   ! that used in the qe2pert_ed calculation
   ftemper = 'si.temper'
/

Notice that an additional variable ftemper specifies the file that contains the doping concentrations and the temperatures, which must be the ones used in the qe2pert_ed interface calculations for charged defects; otherwise, the code would crash. Therefore, one can simply copy the file si.temper used in the qe2pert_ed calculations for charged defects. If one would like to run other different conditions (e.g., temperatures and doping concentations) specified in the si.temper file, one can use the variable lmelonlychar in the qe2pert_ed calculations, and rerun the qe2pert_ed calculation modes “inv_scr”, “edmel” , and “edwan”; note that we do not have to recompute the interpolation table as well as the e-d matrix elements for the neutral part of the defect (i.e., computationally cheap).

calculation mode: ‘setup’

The purpose is to set up the information for transport calculations.

An input file example (“pert.in”) for neutral and charged defects is shown below. The input file is the same for both neutral and charged defects.

&perturbo
  calc_mode = 'setup'
  prefix = 'si'
  tmp_dir = './tmp'

  leph = .false.
  led = .true.

  boltz_kdim(1) = 100
  boltz_kdim(2) = 100
  boltz_kdim(3) = 100

  boltz_emin = 6.0
  boltz_emax = 6.6
  band_min = 2
  band_max = 4
  hole = .true.

  ftemper  = 'si.temper'
/

Here we use a fine k grid of $100\times100\times100$ points. The file ftemper is required for this calculation mode. For neutral defects, the file can contain any doping concentrations and temperatures that users desire to look into; however, for charged defects, the file must be the same one as used in the qe2pert_ed calculations. To use different conditions, please see the description in the section in the above “edmat” calculation mode. The output files will be used for the further transport calculations. The output files are similar to those obtained in the e-ph calculations.

calculation mode: ‘imsigma’

The purpose is to compute the e-d self-energy for the given doping concentration(s) and temperature(s) in the ftemper file. The self-energy can be computed using a regular grid or a random grid, specified by the variable grid_type. Notice that here we turn off leph and turn on led, so the code computes only the e-d self-energy.

An input file example for neutral and charged defects is shown below. Both neutral and charged defects have the same input file. Here we use the regular grid for the final k’ points, and the intial k points are specified in the fklist file. In this case, we use those k points generated in the “setup” calculation modes.

&perturbo
  calc_mode = 'imsigma'
  prefix = 'si'
  tmp_dir = './tmp'

  leph = .false.
  led = .true.

  boltz_kdim(1) = 100
  boltz_kdim(2) = 100
  boltz_kdim(3) = 100

  boltz_emin = 6.0
  boltz_emax = 6.6
  band_min = 2
  band_max = 4
  hole = .true.

  fklist = 'si_tet.kpt'
  ftemper  = 'si.temper'

  ! for a regular grid 
  grid_type = 'regular'

  !! for a random grid 
  !grid_type = 'random'
  !sampling = 'uniform'
  !nsamples = 1000000

  !! for a random grid using the Cauchy sampling used for charged defect
  !grid_type = 'random'
  !sampling = 'cauchy'
  !cauchy_scale = 0.05
  !nsamples = 1000000

  ! If lgeo = .true. (default: .false.), to compute the transport 
  ! scattering rates. 
  !lgeo =.true.
/

To use the uniform regular k’ grid, we need to specify the variable grid_type as “regular”. The output file prefix.imsigma_ed (in this case, si.imsigma_ed) conatins the self-energy for a particular condition (doping and temperature). If the variable lgeo is .true. (the default is .false.), the code also computes the transport scattering rates (which contains the cosine factor), and outputs an additional file called prefix.imsigma_ed_geo (in this case, si.imsigma_ed_geo). Remeber that for charged defects, the file ftemper must be the same as the one used in the interface calculations; to use different conditions, please see the description in the section of the above “edmat” calculation mode.

To use the uniform random k’ grid, we need to specifiy three variables: 1) grid_type as “random”, 2) sampling as “uniform”, and 3) nsamples as the number of final k’ points.

To use the cauchy random k’ grid usually for charged defects, we need to specifiy four variables: 1) “grid_type” as “random”, 2) “sampling” as “cauchy”, 3) “cauchy_scale” as the broadening used in the Cauchy sampling, and 4) “nsamples” as the number of final k’ points.

calculation mode = ‘meanfp’

The purpose is to compute the e-d mean free path using the imsigma file obtained in the calculation mode “imsigma”. For example, in this case, we are inside the directory pert-meanfp, and there we soft link the file si.imsigma_ed using the following command:

# if we use the scattering rates: imsigma_ed
>> ln -sf ../pert-imsigma/si.imsigma_ed ./
# or if we use the transport scattering rates: imsigma_ed_geo
>>  ln -sf ../pert-imsigma/si.imsigma_ed_geo ./si.imsigma_ed
# notice that we need to change the name to imsigma_ed

An input file example for the mean-free-path calculation is as follows:

&perturbo
  prefix = 'si'
  tmp_dir = './tmp'

  leph = .false.
  led = .true.

  calc_mode = 'meanfp'

  boltz_kdim(1) = 100
  boltz_kdim(2) = 100
  boltz_kdim(3) = 100

  boltz_emin = 6.0
  boltz_emax = 6.6
  band_min = 2
  band_max = 4
  hole = .true.

  fklist = 'si_tet.kpt'
  ftemper  = 'si.temper'
/

The format of the output files is the same as that in the e-ph calculation.

calculation mode = ‘trans’

The purpose of this calculation mode is to compute the electrical conductivity or the carrier mobility under the given temperatures and doping concentrations. The input files are basically the same as used in the calculation mode “setup”, except other variables that specifiy the transport calculations such as delta_smear or boltz_nstep. In the following, we introduce three options to compute the transport properties: 1) ITA, 2) RTA, and 3) tr-RTA.

To perform the “trans” calculation mode, one needs to run the “setup” mode first. The “trans” calculation mode requires the following files from the “setup” mode: in this case, si.temper and si_tet.h5, so remember to soft link to those files before running.

ITA appraoch

&perturbo
  calc_mode = 'trans'
  prefix = 'si'
  tmp_dir = './tmp'

  leph = .false.
  led = .true.

  boltz_kdim(1) = 100
  boltz_kdim(2) = 100
  boltz_kdim(3) = 100

  boltz_emin = 6.0
  boltz_emax = 6.6
  band_min = 2
  band_max = 4
  hole = .true.

  ftemper  = 'si.temper'

  delta_smear = 2
  boltz_nstep = 1000

  !to load the calculated e-d matrix elements if they are calculated
  !load_scatter_ed = .true.
/

The description of the output files is the same as that of those in the e-ph calculations. To run the transport properties due to both e-ph and e-d interactions, one just needs to provide the HDF5 files, prefix_epwan.h5 and prefix_edwan.h5 in the directory, and turns on both leph and led.

RTA or tr-RTA approach

This approach requires additional file from the “imsigma” calculation mode: in this case, si.imsigma_ed or si.imsigma_ed_geo. For the later one, make sure that when you soft link to that file, the name of the file shown in the pert-trans directory must be si.imsiga_ed. If the “geo” file is used, the approach is called tr-RTA.

&perturbo
  calc_mode = 'trans'
  prefix = 'si'
  tmp_dir = './tmp'

  leph = .false.
  led = .true.

  boltz_kdim(1) = 100
  boltz_kdim(2) = 100
  boltz_kdim(3) = 100

  boltz_emin = 6.0
  boltz_emax = 6.6
  band_min = 2
  band_max = 4
  hole = .true.

  ftemper  = 'si.temper'

  delta_smear = 2

  ! for RTA and tr-RTA calculations
  ! for RTA, it means we use prefix.imsigma_ed file
  ! for tr-RTA, it means we use prefix.imsigma_ed_geo file 
  boltz_nstep = 0
/

The description of the output files are the same as those in the e-ph calculations.