# GaN-II¶

Author: Samuel Poncé

Warning

The following tutorial is based on an automatic test from the test-suite and is intended to be run quickly. Careful convergence is required to obtain physically meaningful results.

Note

This is a new tutorial which include code features not yet available in the stable version of EPW/QE. Therefore you need the developer version to run it at the moment:

```
git clone https://gitlab.com/QEF/q-e.git
```

## Linearised Boltzmann transport equation (BTE)¶

In this tutorial, we will show how to compute the intrinsic electron and hole mobility of hexagonal wurtzite GaN using the linearised iterative Boltzmann transport equation (IBTE) and the self-energy relaxation time approximation (SERTA).

The theory related to this tutorial can be found in the review First-principles calculations of charge carrier mobility and conductivity in bulk semiconductors and two-dimensional materials.

In particular the IBTE mobility refers to Eqs. 40, 41, 42 and 43 from the review while the SERTA mobility refers to Eqs. 41, 51 and 53. Note that an alternative formulation of the IBTE mobility is also given by Eqs. 63 and 64 by removing the magnetic field (B=0 limit).

In addition, the paper Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors is the recommanded citation for the transport implementation in EPW.

## Preliminary calculations¶

The tutorial files can be found in `q-e/test-suite/epw_mob_polar`

which also include the pseudopotentials.

We first compute the dynamical matrix and electron-phonon matrix element using density functional perturbation theory on the irreducible Brillouin-Zone.

To do so, we first compute the ground-state electron density using 4 cores:

```
mpirun -np 4 ../../bin/pw.x < scf.in | tee scf.out
```

Note

For polar materials as is the case of GaN, the dielectric constant converges very slowly with the k-point grid (linearly). As a result, the k-point grid should be very dense (e.g. 20x20x20). In this example we use 2x2x2 for fast calculation.

We then compute the phonons:

```
mpirun -np 4 ../../bin/ph.x < ph.in | tee ph.out
```

Try looking in the ouput and locate the dielectric constants. You should get 23.12 in-plane and 17.23 for the out of plane direction. Those values are strongly overestimated due to the coarse k-point grid. The calculation should take about 5 min.

Note

The input variable responsible to produce the electron-phonon matrix element is `fildvscf`

. Always make sure that this variable is present.

Finally, we need to post-process some of the data to make it ready for EPW. To do so, we can use the following python script

```
python ../../EPW/bin/pp.py
```

The script will ask you to enter the prefix used for the calculation. In this case enter “gan”. The script will create a new folder called “save” that contains the electron-phonon matrix elements and dynamical matrices on the IBZ.

The script should be compatible with python 2 and python 3.

Note

If you are using python 2, the script might ask you to install the “future” library. If so, type:

```
pip install --user future
```

For information, the converged phonon bandstructure of GaN should look like this:

## Interpolation of the electron-phonon matrix element in real-space with EPW¶

Because the earlier phonon calculation modified the density file, we have to do another DFT self-consistent calculation

```
mpirun -np 4 ../../bin/pw.x < scf.in | tee scf.out
```

Followed by a non-self consistent calculation

```
mpirun -np 4 ../../bin/pw.x -npool 4 < nscf.in | tee nscf.out
```

The reason for the non-self consistent calculation is that EPW needs the wavefunctions on the full BZ on a grid between 0 and 1.

Note

For the nscf calculation the number of pool `-npool`

has to be the same as the total number of core `-np`

Note

Since we are also interested in electron mobility, we will need the conduction bands. Notice that we added the
input `nbnd = 30`

in `nscf.in`

We then run EPW:

```
mpirun -np 4 ../../EPW/src/epw.x -npool 4 < epw1.in | tee epw1.out
```

Note

The number of cores and pool have to be the same as for the `nscf.in`

run.

The following EPW input variables are important in our case:

```
vme = .true. ! Compute exact velocities using Wannier interpolation
lpolar = .true. ! This is a polar material
lifc = .false. ! Do not use external IFC
asr_typ = 'simple' ! Use the simple acoustic sum rule
use_ws = .true. ! Use the optimal Wigner-Seitz cell
nbndsub = 14 ! We Wannierize only 14 bands
nbndskip = 12 ! We skip the first 12 bands
dvscf_dir = './save/' ! Directory where the save folder generated by the @pp.py@ python script is stored
band_plot = .true. ! We want to plot the electronic and/or phonon bandstructure (i.e. we do not use homogeneous fine grids)
filkf = './MGA.txt' ! K-point path from M to Gamma to A
filqf = './MGA.txt' ! Q-point path from M to Gamma to A
efermi_read = .true ! Because our fine grids are non homogeneous, the code cannot determine the Fermi level
fermi_energy= 11.80 ! We need to provide the Fermi level manually.
fsthick = 100.0 ! States that we want to take around the Fermi level. We want all states so we put a large value.@]
```

Note

Because our k and q coarse grid are so coarse, we need `use_ws`

. However it makes the calculation heavier (more vectors).
If your fine grids are big enough (6x6x6 or more), you probably can use Gamma centered WS cell and use `use_ws = .false.`

Note

You can provide the k/q-point path either in crystal or Cartesian coordinate. To do that, modify the second element of the first line of the
`MGA.txt`

file.

Note

To have an idea of which Fermi level you should use, you can first do a calculation with a homogeneous fine grids by removing `filkf`

and `filqf`

and placing:

```
nk1 = 12
nk2 = 12
nk3 = 12
nq1 = 12
nq2 = 12
nq3 = 12
```

Finally, notice that at the end of `epw1.in`

you need to provide the IBZ q-point list:

```
4 cartesian
0.000000000 0.000000000 0.000000000
0.000000000 0.000000000 -0.306767286
0.000000000 -0.577350269 0.000000000
0.000000000 -0.577350269 -0.306767286@]
```

This list has to be exactly the same as the one reported in the output of the phonon calculation at the begining of `ph.out`

that we did earlier.

The calculation should take about 1 min to run. At the end, the code should have produce two files `band.eig`

and `phband.freq`

which contain
the interpolated electronic and phonon bandstructure along the M-G-A path.
You can post-process those file yourself or used the `plotband`

tool from Quantum Espresso by typing

```
../../bin/plotband.x
```

Then follow the instruction. Read the `plotband.x`

documentation if you have an issue. This will create a file which is
graphable with `gnuplot`

and should give you the following:

and the phonon bandstructure:

Given that we have been using a 2x2x2 k-point and q-point grid and a low energy cut-off, such bandstructure quality is remarkable. Comparing with Figure 1, one can see that the only qualitative issue is related to the highest phonon-mode along the G-A direction, which improves with more converged results.

The electronic bandstructure can be compared to the following converged one:

The main difference with the calculation that we just did is that we did not include SOC in this tutorial for speed reasons. You can add SOC by adding

```
noncolin = .true.,
lspinorb = .true.
```

in the `scf.in`

and `nscf.in`

files.

## Calculating the intrinsic drift mobility tensor of GaN with temperature¶

We are now ready to compute the electron and hole mobility tensor with temperature. For this, run the following command to compute the hole mobility:

```
mpirun -np 8 ../../EPW/src/epw.x -npool 8 < epw2.in | tee epw2.out
```

In this case, the code does a restart by reading the electron-phonon matrix element in real-space from the file `gan.epmatwp1`

and interpolate the
matrix element on a fine k/q-point grid.

Note

In case of such a restart, you can use an arbitrary number of cores (.e.g. we used 8 instead of 4). The maximum you can use corresponds to the maximum number of k-points on the fine grid.

Before looking at the output, let us first discuss the specific input variable to perform transport calculations:

```
kmaps = .true. ! Required for a restart
epwread = .true. ! Tels the code to read the gan.epmatwp1 file
wannierize = .false. ! The Wannierization cannot be redone for a restart
scattering = .true. ! We wish to compute scattering rates
iterative_bte = .true. ! We wish to do IBTE and not just SERTA
mob_maxiter = 200 ! Maximum number of iterations for the IBTE
broyden_beta = 1.0 ! Linear mixing between iterations
int_mob = .false. ! We do not want to do electron and hole at the same time
carrier = .true. ! We want to specify carrier concentration
ncarrier = -1E13 ! We want 1E13 carrier per cm^3. The negative sign means hole
mp_mesh_k = .true. ! Use k-point symmetry to reduce the calculation cost
restart = .true. ! We want to create restart point during the interpolation
restart_freq = 50 ! A restart point is created every 50 q-points
selecqread = .false. ! If .true. allows for reading the q-points within the fsthick
epmatkqread = .false. ! Once the transition probabilities have been written to file, allow for a restart to compute the mobilities
nstemp = 2 ! We are considering 2 temperatures.
tempsmin = 100 ! Initial temperature
tempsmax = 300 ! Final temperature
degaussw = 0.00 ! Broadening of the Dirac deltas in the BTE. A value of 0 means adaptative broadening.
```

Note

There is a starting global and unique Fermi level given by `fermi_energy`

or computed by the code.
This Fermi level is used for initialization of quantities and energy window determination.
However the real Fermi level used for the transport calculation is determined automatically by the code depending on the
required carrier concentration and temperature. There can even be up to two such real Fermi level (for electron and hole).

We can take a look at the output of this second EPW calculation. We see the following:

```
Using uniform q-mesh: 12 12 12
Size of q point mesh for interpolation: 1728
Using uniform MP k-mesh: 12 12 12
Size of k point mesh for interpolation: 266
Max number of k points per pool: 34
```

which tells us that we have a q-point grid with 1728 q-points (no symmetry) and a k-point grid with 266 k-points (IBZ). Since we are doing a parallelization on k-point, we will have 34 k-points per core.

Next, the code will try to reduce the number of q-points we need to compute by taking the one that contribute to the
energy window. A file called `selecq.fmt`

will be created and can be read later to avoid recomputing.

```
Number selected, total 50 148
Number selected, total 100 1443
Number selected, total 150 1715
We only need to compute 162 q-points
```

The code then reports the VBM and CBM as well as the Fermi level for the hole calculation for each temperature:

```
Valence band maximum = 11.381872 eV
Conduction band minimum = 13.076046 eV
Temperature 100.000 K
Mobility VB Fermi level = 11.517589 eV
Temperature 300.000 K
Mobility VB Fermi level = 11.804356 eV
```

The code also reports which states are included given the energy window:

```
Fermi Surface thickness = 0.400000 eV
This is computed with respect to the fine Fermi level 11.481872 eV
Only states between 11.081872 eV and 11.881872 eV will be included
```

The code then computes the transition probabilities but only saves the one that are large enough:

```
Save matrix elements larger than threshold: 0.372108862978E-23
```

Note that the threshold is computed automatically based on the size of the system and is very conservative. The total error made by neglecting those matrix elements is lower than machine precision.

Because we are using adaptive smearing, the code reports the min and maximum value of broadening for that q-point.
This information is printed every `restart_freq`

.

```
Progression iq (fine) = 50/ 162
Adaptative smearing = Min: 1.414214 meV
Max: 293.598353 meV
Progression iq (fine) = 100/ 162
Adaptative smearing = Min: 28.908200 meV
Max: 319.312222 meV
Progression iq (fine) = 150/ 162
Adaptative smearing = Min: 48.994246 meV
Max: 149.060946 meV
```

At the end all the transition probabilities are stored in the file `gan.epmatkq1`

as well as the support files `sparsei`

,
`sparsej`

, `sparsek`

, `sparseq`

and `sparset`

. Note that the file size is 28 Kb, orders of magnitude smaller than
without matrix element threshold.

Note

The interpolation and storing of transition probabilities to file can take a lot of time. If the calculation crashes
during runtime (due to time-out for example), just put `selecqread == .true. `` to avoid recomputing the ``selecq.fmt`

and relaunch. The code will automatically continue.

Once all the transition probabilities are compute, you can always restart the calculation by using the `epmatkqread == .true.`

Note that here the code did that for you automatically:

`epmatkqread automatically changed to .TRUE. as all scattering have been computed.`

Then the code computes the SERTA mobility:

```
=============================================================================================
Temp Fermi Hole density Population SR Hole mobility
[K] [eV] [cm^-3] [h per cell] [cm^2/Vs]
=============================================================================================
100.000 11.5176 0.10000E+14 0.12925E-25 0.767956E-02 0.000000E+00 0.000000E+00
0.000000E+00 0.767956E-02 0.000000E+00
0.000000E+00 0.000000E+00 0.573540E+04
300.000 11.8044 0.10000E+14 0.00000E+00 0.585711E+01 0.000000E+00 0.000000E+00
0.000000E+00 0.585711E+01 0.000000E+00
0.000000E+00 0.000000E+00 0.184230E+03
```

The GaN in-plane mobility (xx and yy components) is 0.767956E-02 and the out-of plane one (zz) is 0.573540E+04. Those values are totally not converged and physically wrong because of the very coarse fine grids that we used. Typical grids are larger than 100x100x100. One would always expect the mobility at higher temperature to be lower than at low temperature. This is correctly the case of the out-of plane component but not for the in-plane.

Note

The population sum-rule is simply the sum of all change of population due to external electric field which should be 0. If this value is large, there is probably something wrong in the calculation.

Finally, the code starts solving the IBTE and converges after 36 iterations to:

```
Iteration number: 36
=============================================================================================
Temp Fermi Hole density Population SR Hole mobility
[K] [eV] [cm^-3] [h per cell] [cm^2/Vs]
=============================================================================================
100.000 11.5176 0.10000E+14 -0.12824E-25 0.896091E-02 0.000000E+00 0.000000E+00
0.000000E+00 0.896091E-02 0.000000E+00
0.000000E+00 0.000000E+00 0.374318E+04
300.000 11.8044 0.10000E+14 0.13235E-22 0.757728E+01 0.000000E+00 0.000000E+00
0.000000E+00 0.757728E+01 0.000000E+00
0.000000E+00 0.000000E+00 0.185868E+03
0.806237E-06 Max error
```

To conclude this tutorial, we launch the last calculation with:

```
mpirun -np 8 ../../EPW/src/epw.x -npool 8 < epw3.in | tee epw3.out
```

In this case, we compute the electron mobility.
Take some time to analyse the differences between `epw3.in`

and `epw2.in`

to understand how to do such calculations.
It should be quite self-explanatory.

For information, by using much denser grids as well as SOC and GW corrections on the eigenvalues, one can get a fair agreement with experiment as shown on the following figure:

## Suggestions for efficient transport calculations¶

Transport calculation are computationally heavy as they require a very fine sampling of both the k-point and q-point grids. However, various schemes have been implemented in EPW to reduce the computational cost.

1) The memory scales linearly with the number of temperature `nstemp`

. Try decreasing it if you have memory issues.
The cpu times scales strongly sub-linearly with the number of temperature. Therefore try to do as many temperatures as your memory allows.

2) Low temperatures span a smaller energy window and are therefore more difficult to converge. You can decrease `fsthick`

to make the
calculation more efficient. Note however that a single `fsthick`

is used for all temperature so that your `fsthick`

needs to be large
enough to accommodate the largest temperature.

3) Using `int_mob == .true.`

is convenient because it allows for the calculation of both electron and hole at the same time.
However, the code consider only the element within an energy window given by `fsthick`

around the starting Fermi level given by `fermi_energy`

. As a result you might need to use a larger than needed `fsthick`

. The recommended way is therefore to compute
electron and hole mobility in two separate runs.

4) Transport properties only require states very close to the band edge and significant speedup can be achieved
by using a small `fsthick`

. Note that you need to converge your result with increasing values of `fsthick`

and that larger values
are required for higher temperature. For example in this tutorial we used a 0.4 eV `fsthick`

but we placed our Fermi level
0.1 eV above the valence band maximum (VBM). As a result we are considering only the states up to 0.3 eV above the VBM.
This is a typical order of magnitude.