Part 5: Phonon-limited mobility of semiconductors and the ZT figure of merit ¶
By the end of this tutorial, you will be able to:
- Compute mobilities of silicon in the constant relaxation-time approximation (CRTA).
- Perform numeric convergence tests of energy integration grids and $\textbf{k}$-points to calculate the electron mobility.
- Compare the two available approaches to integrate the transport distribution function and obtain the Onsager coefficients.
14.1 Task¶
Calculate the electron mobility for cd Si in the constant relaxation-time approximation (CRTA) and converge for the energy grid (linear and Gauss-Legendre) used to calculate the Onsager coefficients, and k-point mesh.
Unlike for metals, where many charge carriers are present due to the zero bandgap, semiconductors have no charge carriers at $T \neq 0$. Carriers are added by doping the semiconductor to increase the charge density (electrons or holes for n- and p-type, respectively). The ease with which the charge carriers are able to move through the material is described by the mobility $\mu_e$ or $\mu_h$ (electron or hole), a key indicator of the performance of semiconductors in electrical devices, related to their conductivity $\sigma$, which can be defined using the transport distribution function $\mathcal{T}(\epsilon)$ (cf. Eq. 11.1).
By restricting the band sums to conduction or valence bands, the conductivity can be separated into electron and hole contributions, $\sigma_e$ and $\sigma_h$, and the electron and hole mobilities defined:
| Quantity | Definition | Carrier density |
|---|---|---|
| Electron mobility $\mu_e$ | $\mu_e = \tfrac{\sigma_{n \in \text{CB}}}{n_e}$ | $n_e = \frac{1}{\Omega N_\mathbf{k}} \sum_{\mathbf{k}n \in \text{CB}} f(\varepsilon, T, \mu)$ |
| Hole mobility $\mu_h$ | $\mu_h = \tfrac{\sigma_{n \in \text{VB}}}{n_h}$ | $n_h = \frac{1}{\Omega N_\mathbf{k}} \sum_{\mathbf{k}n \in \text{VB}} \big[1 - f(\varepsilon, T, \mu)\big]$ |
where $\sigma_{n \in \text{CB}}$ and $\sigma_{n \in \text{VB}}$ denote the conductivity restricted to states in the conduction and valence bands, respectively, $f_{n\mathbf{k}}$ is the Fermi–Dirac distribution, and $\mu$ is the chemical potential at the given temperature.
Note, the chemical potential is always written as $\mu$. The mobility always has a subscripts, e.g.: $\mu_e$ or $\mu_h$.
In this tutorial, you will explore how to compute the mobility of electrons in semiconductors using the linearized Boltzmann transport equation (BTE), with a focus on the effects of electron-phonon interactions [Rep. Prog. Phys. 83 036501 (2020)]. In the CRTA, the mobility is only impacted by the geometry of the electronic band structure (i.e., the carriers heaviness/ lightness, the anisotropy of the valleys, and reordering due to strain). The CRTA mobility is a “band-structure fingerprint”, indicating how the crystal would behave without scattering.
14.2 Transport coefficients in the constant relaxation-time approximation¶
14.2.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e14_mobility_si_crta/.
Silicon
1.0000000000000000
2.7149999999999999 -2.7149999999999999 0.0000000000000000
2.7149999999999999 0.0000000000000000 2.7149999999999994
0.0000000000000000 -2.7149999999999999 2.7150000000000003
Si
2
Direct
0.2500000000000000 0.2500000000000000 0.2500000000000002
0.5000000000000001 0.4999999999999999 0.4999999999999999
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
TRANSPORT_RELAXATION_TIME=1e14
Automatic
0
Gamma
8 8 8
Automatic
0
Gamma
24 24 24
PAW_PBE Si 05Jan2001
The POSCAR file defines the primitive cell of silicon.
The ELPH_SCATTERING_APPROX=CRTA tag means that we do not need to compute the electron-phonon potential in a supercell beforehand. CRTA is very crude for semiconductors, as you will see in the next section, because of the sharp behaviour of the density of states. The tag ELPH_TRANSPORT_DRIVER=2 specifies to use Gauss-Legendre integration for the energy integration of the transport distribution function. ELPH_SELFEN_TEMPS gives the temperatures at which to compute the self-energy due to electron-phonon scattering. ELPH_ISMEAR=-24 chooses the smearing method to determine the Fermi level and chemical potential before an electron-phonon calculation. -24 is a more memory- and time-efficient implementation of the tetrahedron method (you can try setting ELPH_ISMEAR=-15, and it will take ~6-7 times longer). ELPH_SELFEN_CARRIER_DEN lists the number of additional carriers. Note: ELPH_SELFEN_CARRIER_DEN adds electrons when the numbers are positive n-doping (to add electrons), and decreases when negative p-doping (to add holes), the reverse of the literature convention.
14.2.2 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e14_*/crta
mpirun -np 4 vasp_std
This will run the same calculation of the eigenvalues and electron group velocities and evaluate the transport distribution function for different Gauss-Legendre energy grids. It is important to converge the transport coefficients with these energy grids, in particular, their interplay with ELPH_SELFEN_CARRIER_DENSITY.
At this point, the cost of different grid densities is small because we use the CRTA. Once the calculation is finished, we can plot the mobility as a function of carrier density for the different energy grid samplings we choose.
import pandas as pd
import py4vasp
import plotly.graph_objects as go
from py4vasp import Calculation
calculation = Calculation.from_path('e14_mobility_si_crta/crta')
mobility = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][:,0,0]
temperatures = transport.temperatures()
mobility[carrier_den] = dict(zip(temperatures, mob))
df = pd.DataFrame(mobility)
df = df.sort_index()
fig = go.Figure()
for temp in df.index:
fig.add_trace(go.Scatter(
x=df.columns,
y=df.loc[temp],
mode='markers+lines',
name=f'{temp} K',
marker=dict(size=8),
line=dict(width=2)
))
# Layout
fig.update_layout(
title='Electron mobility in Si for different temperatures',
xaxis=dict(
title='Carrier density in cm<sup>-3</sup>',
type='log',
showgrid=True,
gridwidth=1
),
yaxis=dict(
title='Mobility in cm<sup>2</sup>/Vs',
showgrid=True,
gridwidth=1
),
width=800,
height=500
)
fig.show()
The mobility is dependent on the carrier density. At high densities, this is correct physical behaviour; at low densities, the mobility should be independent of the density. This hump is indicative that convergence should be done with respect to ELPH_TRANSPORT_NEDOS.
The value of the mobility using CRTA is poor when compared with experiment. While the experimental values are around 1400 cm/Vs for electrons and 500 cm/Vs for holes [Phys. Rev. B 97 121201(R) (2018)], using CRTA we get about 99 cm/Vs and 72 cm/Vs. This is to be expected, since we have set the relaxation time at a constant 1E-14. Tuning this parameter would allow us to recover the correct experimental value.
14.3 Linear grids and Simpson integration¶
The Onsager coefficients $L_{ij}$ can be used to calculate the transport coefficients. They are defined as: $$\tag{14.1} L_{ij} = \int d\epsilon \, \mathcal{T}(\epsilon) \, (\epsilon-\mu)^{i+j-2} \left( -\frac{\partial f^0}{\partial \epsilon} \right), $$ where $\mathcal{T}(\epsilon)$ is the transport distribution function, $\mu$ the chemical potential, and $f^0$ the Fermi–Dirac distribution, and are integrated using Simpson's integration.
In this section, we will deal with details of how the energy integral of the transport distribution function is performed numerically and the transport coefficients obtained.
14.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e14_mobility_si_crta/linear_grid.
51/INCAR
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_TRANSPORT_DRIVER=1
ELPH_TRANSPORT_DFERMI_TOL=1e-6
ELPH_FERMI_NEDOS=51
ELPH_TRANSPORT_NEDOS=51
Click to see 501/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_TRANSPORT_DRIVER=1
ELPH_TRANSPORT_DFERMI_TOL=1e-6
ELPH_FERMI_NEDOS=501
ELPH_TRANSPORT_NEDOS=501
Click to see 5001/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_TRANSPORT_DRIVER=1
ELPH_TRANSPORT_DFERMI_TOL=1e-6
ELPH_FERMI_NEDOS=5001
ELPH_TRANSPORT_NEDOS=5001
Click to see 10001/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_TRANSPORT_DRIVER=1
ELPH_TRANSPORT_DFERMI_TOL=1e-6
ELPH_FERMI_NEDOS=10001
ELPH_TRANSPORT_NEDOS=10001
Automatic
0
Gamma
8 8 8
Automatically generated mesh
0
Gamma
24 24 24
PAW_PBE Si 05Jan2001
There are only a few unfamiliar INCAR tags. ELPH_TRANSPORT_DRIVER=1 selects the linear energy grid to evaluate the Onsager coefficients, ELPH_FERMI_NEDOS and ELPH_TRANSPORT_NEDOS set the number of linear energy grid points, and ELPH_TRANSPORT_DFERMI_TOL defines the energy window for calculating the Onsager coefficients.
14.3.2 Calculation¶
This will run the same calculation of the eigenvalues and electron group velocities and evaluate the transport distribution function for different numbers of points in the energy grids. The goal is to monitor the impact that these energy grids can have on the numerical accuracy of the transport coefficients, and in particular, their interplay with ELPH_SELFEN_CARRIER_DENSITY.
We used the CRTA to test the numeric differences between grids, as the calculations are quick. A converged grid for CRTA is the smallest you should try for SERTA.
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e14_*/linear_grid
bash linear_grid.sh
Click to see linear_grid.sh!
for nedos in 51 501 5001 10001
do
cd $nedos
mpirun -np 4 vasp_std
cd ..
done
Once the calculation is finished we can plot the mobility as a function of carrier density for the different energy grid samplings we choose:
import pandas as pd
import plotly.graph_objects as go
from py4vasp import Calculation
TRANSPORT_NEDOS = [51, 501, 5001, 10001]
itemp = 3
mobility_kpoints = {}
for n in TRANSPORT_NEDOS:
calculation = Calculation.from_path(f'e14_mobility_si_crta/linear_grid/{n}')
mobility = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][itemp][0][0]
if carrier_den < 0:
continue
mobility[carrier_den] = mob
mobility_kpoints[f'{n}'] = mobility
df = pd.DataFrame(mobility_kpoints)
df = df.sort_index()
fig = go.Figure()
for col in df.columns:
fig.add_trace(go.Scatter(
x=df.index,
y=df[col],
mode='markers+lines',
name=f'TRANSPORT_NEDOS = {col}',
marker=dict(size=8),
line=dict(width=2)
))
fig.update_layout(
title='Electron mobility in CRTA using linear grids for Si at 300 K',
xaxis=dict(
title='Carrier density in cm<sup>-3</sup>',
type='log',
showgrid=True,
gridwidth=1
),
yaxis=dict(
title='Mobility in cm<sup>2</sup>/Vs',
showgrid=True,
gridwidth=1
),
width=800,
height=500
)
fig.show()
You can see how the blue, cyan, and red lines coincide, indicating energy grid convergence has been reached for each carrier density.
The main drawback of this approach is that the energy range is chosen by the ELPH_TRANSPORT_DFERMI_TOL parameter. Specifically, the energy window is chosen so that a certain fraction of the integral of the derivative of the Fermi-Dirac distribution around the chemical potential is excluded. As a result, simply increasing the number of energy points is not sufficient to ensure convergence of the mobility; it is also necessary to converge with respect to ELPH_TRANSPORT_DFERMI_TOL.
As an optional exercise, try changing ELPH_TRANSPORT_DFERMI_TOL=1e-6 to ELPH_TRANSPORT_DFERMI_TOL=1e-8 in the INCAR files above, then re-run the calculation and plot the results to observe the effect.
With a stricter tolerance, less of the integral is neglected, and so the mobility should not depend on the carrier density at low densities. The bump before 1016 indicates an issue with numerical convergence.
To avoid this inconvenient double convergence, the Gauss-Legendre integration method has been implemented for the transport coefficients. In the next section, we will use Gauss-Legendre grids to compute the mobilities and compare the results.
14.4 Gauss-Legendre integration¶
14.4.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e14_mobility_si_crta/gauss_legendre.
51/INCAR
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_FERMI_NEDOS=51
ELPH_TRANSPORT_NEDOS=51
Click to see 501/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_FERMI_NEDOS=501
ELPH_TRANSPORT_NEDOS=501
Click to see 5001/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_FERMI_NEDOS=5001
ELPH_TRANSPORT_NEDOS=5001
Click to see 10001/INCAR!
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_FERMI_NEDOS=10001
ELPH_TRANSPORT_NEDOS=10001
Automatic
0
Gamma
8 8 8
Automatically generated mesh
0
Gamma
24 24 24
PAW_PBE Si 05Jan2001
Now we compute again the mobility in the CRTA approximation but this time using Gauss-Legendre integration. This is selected by ELPH_TRANSPORT_DRIVER=2, which is the default.
14.4.2 Calculation¶
Now we will perform the same calculation but using Gauss-Legendre grids.
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e14_*/gauss_legendre
bash gauss_legendre.sh
Click to see gauss_legendre.sh!
for nedos in 51 501 5001 10001
do
cd $nedos
mpirun -np 4 vasp_std
cd ..
done
While the calculations are running, have a look at the documentation of ELPH_TRANSPORT_DRIVER=2 and get yourself acquainted with how the energy grids are chosen and how this changes the energy range at which the transport distribution function is to be evaluated.
You can find out which grid has been used in the electron-phonon section of the OUTCAR file, written at the end of the line for the electron-phonon accumulators as Linear grids:
Transport for self-energy accumulator N= 4
T K mu eV sigma S/m mob cm^2/(V.s) seebeck μV/K peltier μV kappa_e W/(m.K)
0.00000000 5.82498586 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 Linear grids
or Gauss-Legendre grids:
Transport for self-energy accumulator N= 4
T K mu eV sigma S/m mob cm^2/(V.s) seebeck μV/K peltier μV kappa_e W/(m.K)
0.00000000 5.82498587 115.92371848 72.35410320 0.00000000 0.00000000 0.00000000 Gauss-Legendre grid
Once the calculations are finalized, you can plot the electron mobility as a function of carrier doping for the different grids with the following script:
import pandas as pd
import py4vasp
import plotly.graph_objects as go
from py4vasp import Calculation
TRANSPORT_NEDOS = [51, 501, 5001, 10001]
itemp = 3
mobility_kpoints = {}
for n in TRANSPORT_NEDOS:
calculation = Calculation.from_path(f'e14_mobility_si_crta/gauss_legendre/{n}')
mobility = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][itemp][0][0]
mobility[carrier_den] = mob
mobility_kpoints[f'{n}'] = mobility
df = pd.DataFrame(mobility_kpoints)
df = df.sort_index()
fig = go.Figure()
for col in df.columns:
fig.add_trace(go.Scatter(
x=df.index,
y=df[col],
mode='markers+lines',
name=f'TRANSPORT_NEDOS = {col}',
marker=dict(size=8)
))
# Layout
fig.update_layout(
title='Electron mobility in CRTA using Gauss-Legendre grids for Si at 300 K',
xaxis_title='Carrier density in cm<sup>-3</sup>',
yaxis_title='Mobility in cm<sup>2</sup>/Vs',
xaxis=dict(type='log', showgrid=True, gridwidth=1),
yaxis=dict(showgrid=True, gridwidth=1),
width=1000,
height=500
)
fig.show()
You can see that, though the convergence is slower for the number of energy points, there is no need to tune an additional parameter to choose the energy range (ELPH_TRANSPORT_DFERMI_TOL in the case of ELPH_TRANSPORT_DRIVER=1). Reducing the number of parameters from two to one makes it simpler to converge. The linear grid can speed up calculations, but requires more experience and expertise.
This is particularly important when computing the electronic lifetimes from first principles, i.e., when setting, for example, ELPH_SCATTERING_APPROX=SERTA. The grid should be converged using CRTA, as it is the lower limit for SERTA. The computation of the lifetimes is rather time-consuming, as computing the lifetime for each $\tau_{n\mathbf{k}}$ implies a sum over $\mathbf{q}$. We can greatly reduce the amount of computation by computing a smaller subset of $\tau_{n\mathbf{k}}$ for fewer bands and $\textbf{k}$-points. The energy range is used to select which $\tau_{n\mathbf{k}}$ need to be computed, thus greatly reducing the amount of computational effort. The Gauss–Legendre approach has the advantage that the integration grid adapts naturally to the width of the Fermi window, making it numerically efficient without having to manually define the energy window (cf. ELPH_TRANSPORT_DFERMI_TOL, and ELPH_TRANSPORT_EMIN and ELPH_TRANSPORT_EMAX). Instead, only the number of points in the sum above needs to be defined through TRANSPORT_NEDOS.
14.5 Convergence with k-points¶
Now that we have compared then different energy integration methods to go from the transport distribution function to the transport coefficients you are ready to perform a convergence test with $\textbf{k}$-point sampling.
14.5.1 Input¶
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e16 1e16
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=CRTA
#nedos
ELPH_FERMI_NEDOS=5001
ELPH_TRANSPORT_NEDOS=5001
Automatic
0
Gamma
8 8 8
Click to see q-mesh-12x12x12/KPOINTS_ELPH!
Automatic
0
Gamma
12 12 12
Click to see q-mesh-18x18x18/KPOINTS_ELPH!
Automatic
0
Gamma
18 18 18
Click to see q-mesh-24x24x24/KPOINTS_ELPH!
Automatic
0
Gamma
24 24 24
Click to see q-mesh-32x32x32/KPOINTS_ELPH!
Automatic
0
Gamma
32 32 32
Click to see q-mesh-48x48x48/KPOINTS_ELPH!
Automatic
0
Gamma
48 48 48
PAW_PBE Si 05Jan2001
14.5.2 Calculation¶
To run these calculations for the different $\textbf{k}$-point meshes you can copy and paste the following code into a terminal:
cd $TUTORIALS/electron-phonon/e14_*/kpoints
bash kpoints.sh
Click to see kpoints.sh!
for k in 12 18 24 32
do
cd ${k}x${k}x${k}
mpirun -np 4 vasp_std
cd ..
done
Plot the convergence of the electron and hole mobility with respect to $\textbf{k}$-points using the following Python script:
import pandas as pd
import plotly.graph_objects as go
from py4vasp import Calculation
itemp = 3
mobility_kpoints = {}
kmeshes = []
for k in [12, 18, 24, 32, 48, 56]:
calculation = Calculation.from_path(f'e14_mobility_si_crta/kpoints/{k}x{k}x{k}')
mobility = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][itemp][0][0]
if abs(carrier_den - 1e16) < 1:
mobility['electrons'] = mob
if abs(carrier_den + 1e16) < 1:
mobility['holes'] = mob
kmeshes.append(k)
mobility_kpoints[k] = mobility
df = pd.DataFrame(mobility_kpoints).T
fig = go.Figure()
# Holes
fig.add_trace(go.Scatter(
x=df.index, y=df['holes'],
mode='markers+lines',
name='Holes',
marker=dict(symbol='circle', size=8)
))
# Electrons
fig.add_trace(go.Scatter(
x=df.index, y=df['electrons'],
mode='markers+lines',
name='Electrons',
marker=dict(symbol='circle', size=8)
))
# Layout
fig.update_layout(
title='Electron mobility in CRTA for Si at 300 K',
xaxis_title='k-point sampling',
yaxis_title='Mobility in cm<sup>2</sup>/Vs',
xaxis=dict(range=[0, max(df.index)*1.05]),
width=900,
height=600
)
fig.show()
You can see that from an initial $\textbf{k}$-point mesh of 12x12x12, the mobilities of both the holes and the electrons begin to converge, though convergence requires a much denser mesh than we do here.
A more convenient way to reach an "infinite" $\textbf{k}$-point mesh is by plotting the mobility as a function of the inverse of the number of divisions used to generate your $\textbf{k}$-point mesh and extrapolating. In the following script, you can linearly extrapolate the results to infinite $\textbf{k}$-point sampling:
import numpy as np
import pandas as pd
import py4vasp
import plotly.graph_objects as go
df = pd.DataFrame(mobility_kpoints)
inv_k = 1 / np.array(kmeshes)
df_T = df.T
df_T.index = inv_k
# Extrapolations
# Holes: last 3 points
x_hole = inv_k[3:]
y_hole = df_T['holes'].values[3:]
hole_coeffs = np.polyfit(x_hole, y_hole, 1)
hole_fit = np.poly1d(hole_coeffs)
# Electrons: last 3 points
x_elec = inv_k[3:]
y_elec = df_T['electrons'].values[3:]
elec_coeffs = np.polyfit(x_elec, y_elec, 1)
elec_fit = np.poly1d(elec_coeffs)
# Extrapolation x values
x_extrap = np.linspace(0, inv_k.max(), 100)
# Plotly
fig = go.Figure()
# Original data
fig.add_trace(go.Scatter(
x=inv_k, y=df_T['holes'],
mode='markers+lines',
name='Holes',
marker=dict(symbol='circle', size=8)
))
fig.add_trace(go.Scatter(
x=inv_k, y=df_T['electrons'],
mode='markers+lines',
name='Electrons',
marker=dict(symbol='circle', size=8)
))
# Extrapolation lines
fig.add_trace(go.Scatter(
x=x_extrap, y=hole_fit(x_extrap),
mode='lines',
line=dict(dash='dash', color='#89ad01'),
name='Hole extrapolation (3 pts)'
))
fig.add_trace(go.Scatter(
x=x_extrap, y=elec_fit(x_extrap),
mode='lines',
line=dict(dash='dash', color='#A82C35'),
name='Electron extrapolation (3 pts)'
))
# Layout
fig.update_layout(
title='Electron mobility in CRTA for Si at 300 K',
xaxis_title='1 / k-point sampling',
yaxis_title='Mobility in cm<sup>2</sup>/Vs',
width=1000,
height=600
)
fig.show()
From the inverted $\textbf{k}$-point sampling, it is clear that the mobility can be extrapolated to the infinite $\textbf{k}$-mesh.
Congratulations on obtaining converged results! We have been using an ideal, perfect crystal. In a real system, there are also ionized impurities (defects), which affect the scattering rate of the carriers. At high defect density, the mobility of the semiconductor is reduced. In doping, you increase the number of carriers and defects, increasing the carrier density. However, after a certain point, this reaps diminishing returns as the increased number of defects limits the mobility of the charge carriers.
14.5 Questions¶
- What INCAR tags can you use to define the Gauss-Legendre energy grids?
- How can one extrapolate the mobility to infinite $\textbf{k}$-point sampling?
By the end of this tutorial, you will be able to do the following:
- Compute the electron-phonon potential on a supercell.
- Converge the mobility in SERTA with respect to the number of $\textbf{k}$- and $\textbf{q}$-points together.
- Observe how the carrier mobility for holes and electrons changes with $\textbf{k}$ convergence.
15.1 Task¶
Calculate the electron-phonon potential for a 2x2x2 Si supercell, then converge the mobilities with respect to k- and q-point mesh for a cd Si cell in the self-energy relaxation-time approximation (SERTA)
In CRTA, the electron lifetime (relaxation time) is constant. By changing this in SERTA, the mobility will be impacted, due to the inclusion of electron-phonon coupling from first-principles. The first step for this is preparing the electron-phonon potential in a supercell, as you have done previously, followed by a convergence study.
15.2 Electron-phonon potential from a supercell¶
15.2.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e15_mobility_si_serta/supercell.
Each subdirectory specifies a different $\textbf{q}$-point mesh via the KPOINTS_ELPH file.
Check them out!
Click to see POSCAR for the 2x2x2 Si supercell!
2x2x2 Si supercell
1.0
10.8599999999999994 0.0000000000000000 0.0000000000000000
0.0000000000000000 10.8599999999999994 0.0000000000000000
0.0000000000000000 0.0000000000000000 10.8599999999999994
Si
64
Direct
0.1250000000000000 0.1250000000000000 0.3750000000000000
0.6250000000000000 0.1250000000000000 0.3750000000000000
0.1250000000000000 0.6250000000000000 0.3750000000000000
0.6250000000000000 0.6250000000000000 0.3750000000000000
0.1250000000000000 0.1250000000000000 0.8750000000000000
0.6250000000000000 0.1250000000000000 0.8750000000000000
0.1250000000000000 0.6250000000000000 0.8750000000000000
0.6250000000000000 0.6250000000000000 0.8750000000000000
0.2500000000000000 0.0000000000000000 0.0000000000000000
0.7500000000000000 0.0000000000000000 0.0000000000000000
0.2500000000000000 0.5000000000000000 0.0000000000000000
0.7500000000000000 0.5000000000000000 0.0000000000000000
0.2500000000000000 0.0000000000000000 0.5000000000000000
0.7500000000000000 0.0000000000000000 0.5000000000000000
0.2500000000000000 0.5000000000000000 0.5000000000000000
0.7500000000000000 0.5000000000000000 0.5000000000000000
0.1250000000000000 0.3750000000000000 0.1250000000000000
0.6250000000000000 0.3750000000000000 0.1250000000000000
0.1250000000000000 0.8750000000000000 0.1250000000000000
0.6250000000000000 0.8750000000000000 0.1250000000000000
0.1250000000000000 0.3750000000000000 0.6250000000000000
0.6250000000000000 0.3750000000000000 0.6250000000000000
0.1250000000000000 0.8750000000000000 0.6250000000000000
0.6250000000000000 0.8750000000000000 0.6250000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
0.7500000000000000 0.2500000000000000 0.2500000000000000
0.2500000000000000 0.7500000000000000 0.2500000000000000
0.7500000000000000 0.7500000000000000 0.2500000000000000
0.2500000000000000 0.2500000000000000 0.7500000000000000
0.7500000000000000 0.2500000000000000 0.7500000000000000
0.2500000000000000 0.7500000000000000 0.7500000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
0.3750000000000000 0.1250000000000000 0.1250000000000000
0.8750000000000000 0.1250000000000000 0.1250000000000000
0.3750000000000000 0.6250000000000000 0.1250000000000000
0.8750000000000000 0.6250000000000000 0.1250000000000000
0.3750000000000000 0.1250000000000000 0.6250000000000000
0.8750000000000000 0.1250000000000000 0.6250000000000000
0.3750000000000000 0.6250000000000000 0.6250000000000000
0.8750000000000000 0.6250000000000000 0.6250000000000000
0.0000000000000000 0.0000000000000000 0.2500000000000000
0.5000000000000000 0.0000000000000000 0.2500000000000000
0.0000000000000000 0.5000000000000000 0.2500000000000000
0.5000000000000000 0.5000000000000000 0.2500000000000000
0.0000000000000000 0.0000000000000000 0.7500000000000000
0.5000000000000000 0.0000000000000000 0.7500000000000000
0.0000000000000000 0.5000000000000000 0.7500000000000000
0.5000000000000000 0.5000000000000000 0.7500000000000000
0.3750000000000000 0.3750000000000000 0.3750000000000000
0.8750000000000000 0.3750000000000000 0.3750000000000000
0.3750000000000000 0.8750000000000000 0.3750000000000000
0.8750000000000000 0.8750000000000000 0.3750000000000000
0.3750000000000000 0.3750000000000000 0.8750000000000000
0.8750000000000000 0.3750000000000000 0.8750000000000000
0.3750000000000000 0.8750000000000000 0.8750000000000000
0.8750000000000000 0.8750000000000000 0.8750000000000000
0.0000000000000000 0.2500000000000000 0.0000000000000000
0.5000000000000000 0.2500000000000000 0.0000000000000000
0.0000000000000000 0.7500000000000000 0.0000000000000000
0.5000000000000000 0.7500000000000000 0.0000000000000000
0.0000000000000000 0.2500000000000000 0.5000000000000000
0.5000000000000000 0.2500000000000000 0.5000000000000000
0.0000000000000000 0.7500000000000000 0.5000000000000000
0.5000000000000000 0.7500000000000000 0.5000000000000000
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#supercell
IBRION=6
ELPH_POT_GENERATE=TRUE
Automatic
0
Gamma
1 1 1
PAW_PBE Si 05Jan2001
15.2.2 Calculation¶
Since the supercell calculation takes some time (~5 minutes), we will omit this from the exercise. All of the files can be found in the directory, so you can run the calculation if you have time. We provide the phelel_params.hdf5 file for the following electron-phonon calculations.
Click here to see how to run the calculation!
cd $TUTORIALS/electron-phonon/e15_*/supercell
mpirun -np 4 vasp_std
15.3 Convergence with respect to the k- and q-point mesh¶
Now that we have computed the electron-phonon potential on a supercell, we are ready to perform a full first-principles calculation of the mobility. By explicitly choosing ELPH_SCATTERING_APPROX=SERTA, the electron-phonon potential will be read from phelel_params.hdf5 and interpolated onto the $\textbf{q}$-mesh that connects all the $\textbf{k}$-points defined in the KPOINTS_ELPH file. As a result, changing the mesh in the KPOINTS_ELPH file simultaneously adjusts both the $\textbf{k}$- and $\textbf{q}$-point meshes. This allows us to systematically study the convergence of the calculated mobility with respect to the density of these meshes.
15.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e15_mobility_si_serta/kpoints.
Each subdirectory specifies a different $\textbf{q}$-point mesh via the KPOINTS_ELPH file.
Check them out!
Silicon
1.0000000000000000
2.7149999999999999 -2.7149999999999999 0.0000000000000000
2.7149999999999999 0.0000000000000000 2.7149999999999994
0.0000000000000000 -2.7149999999999999 2.7150000000000003
Si
2
Direct
0.2500000000000000 0.2500000000000000 0.2500000000000002
0.5000000000000001 0.4999999999999999 0.4999999999999999
EDIFF=1e-8
ISMEAR=0
SIGMA=0.05
PREC=Accurate
#paralelization
KPAR=4
#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=-1e20 -1e19 -1e18 -1e17 -1e16 -1e15 -1e14 -1e13 1e13 1e14 1e15 1e16 1e17 1e18 1e19 1e20
ELPH_SELFEN_TEMPS=0 100 200 300
ELPH_ISMEAR=-24
ELPH_SCATTERING_APPROX=SERTA
# nedos
ELPH_FERMI_NEDOS=5001
ELPH_TRANSPORT_NEDOS=5001
Automatic
0
Gamma
8 8 8
PAW_PBE Si 05Jan2001
Click to see q-mesh-12x12x12/KPOINTS_ELPH!
Automatic
0
Gamma
12 12 12
Click to see q-mesh-18x18x18/KPOINTS_ELPH!
Automatic
0
Gamma
18 18 18
Click to see q-mesh-24x24x24/KPOINTS_ELPH!
Automatic
0
Gamma
24 24 24
Click to see q-mesh-32x32x32/KPOINTS_ELPH!
Automatic
0
Gamma
32 32 32
Click to see q-mesh-48x48x48/KPOINTS_ELPH!
Automatic
0
Gamma
48 48 48
15.3.2 Calculation¶
Open a terminal and navigate to the kpoints directory of this exercise. Next, run an electron-phonon calculation for each $\textbf{q}$-point mesh by looping over the directories using the qpoints.sh script:
cd $TUTORIALS/electron-phonon/e15_*/kpoints
bash qpoints.sh
Click to see qpoints.sh!
for q in 12 18; do
cp ../supercell/phelel_params.hdf5 "${q}x${q}x${q}"
cd "${q}x${q}x${q}"
mpirun -np 4 vasp_std
cd ../
done
Due to time constraints, we perform only on the 12x12x12 and 18x18x18 grids. The remainder take 10+ minutes so are pre-calculated.
While the calculation is running it is a good opportunity to get yourself aquainted with the different possible approximations for the scattering rate, ELPH_SCATTERING_APPROX.
When it is finished, try plotting the electron mobility against the carrier density:
from py4vasp import Calculation
import py4vasp
import plotly.graph_objects as go
# --- Step 1: Define k-point meshes and temperature index ---
kmeshes = [12, 18, 24, 32, 48]
itemp = 3 # index for 300 K
# --- Step 2: Collect mobility data ---
# Structure: mobility_kpoints[kmesh] = dict of carrier_density -> mobility
mobility_kpoints = {}
for k in kmeshes:
calculation = Calculation.from_path(f'e15_mobility_si_serta/kpoints/{k}x{k}x{k}')
mobility_kpoints[f'{k}x{k}x{k}'] = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][itemp][0][0]
# Skip invalid or negative densities
if carrier_den < 0:
continue
# Store as numeric keys
mobility_kpoints[f'{k}x{k}x{k}'][float(carrier_den)] = float(mob)
# --- Step 3: Plot using Plotly ---
fig = go.Figure()
for i, (kmesh, data_dict) in enumerate(mobility_kpoints.items()):
if not data_dict:
continue # skip empty
# Sort by carrier density
sorted_pairs = sorted(data_dict.items())
x_vals, y_vals = zip(*sorted_pairs)
fig.add_trace(
go.Scatter(
x=x_vals,
y=y_vals,
mode='lines+markers',
name=f"{kmesh}",
marker=dict(size=6),
hovertemplate=(
f"k-mesh: {kmesh}<br>"
"Density: %{x:.2e} cm⁻³<br>"
"Mobility: %{y:.2f} cm²/Vs"
)
)
)
# --- Step 4: Layout and axes ---
fig.update_layout(
title="Electron mobility in Si at 300 K",
xaxis=dict(
title="Carrier density in cm<sup>-3</sup>",
type="log",
showgrid=True
),
yaxis=dict(
title="Mobility in cm<sup>2</sup>/Vs",
showgrid=True
),
legend_title="k-points",
width=900,
height=500
)
fig.show()
Compared to the curve in Section 14.3, you can see that the numerical difficulties have been resolved and the mobility is independent of the carrier density at low densities.
You can see that as the $\textbf{k}$-mesh density increases, the mobility of the electrons also increases until convergence is reached. You can plot the mobility of the holes and electrons separately using the following script:
from py4vasp import Calculation
import py4vasp
import plotly.graph_objects as go
meshes = [12, 18, 24, 32, 48]
itemp = 3
mobility_kpoints = {}
for k in kmeshes:
calculation = Calculation.from_path(f'e15_mobility_si_serta/kpoints/{k}x{k}x{k}')
mobility = {}
for transport in calculation.electron_phonon.transport:
transport_dict = transport.to_dict()
carrier_den = transport_dict['metadata']['selfen_carrier_den']
mob = transport_dict['mobility'][itemp][0][0]
if abs(carrier_den-1e16) < 1: # electrons
mobility['electrons'] = mob
elif abs(carrier_den+1e16) < 1: # holes
mobility['holes'] = mob
mobility_kpoints[k] = mobility
# --- Plotting ---
fig = go.Figure()
for carrier in ['electrons', 'holes']:
x_vals = []
y_vals = []
for k in sorted(mobility_kpoints.keys()):
if carrier in mobility_kpoints[k]:
x_vals.append(k)
y_vals.append(mobility_kpoints[k][carrier])
fig.add_trace(
go.Scatter(
x=x_vals,
y=y_vals,
mode='lines+markers',
name=carrier.capitalize(),
marker=dict(size=8)
)
)
fig.update_layout(
title="Carrier mobility in Si at 300 K",
xaxis=dict(title="k- and q-point sampling (identical grids)"),
yaxis=dict(title="Mobility in cm<sup>2</sup>/Vs"),
legend_title="Carrier",
width=800,
height=500
)
fig.show()
Congratulations on completing another convergence study!
It is clear that the mobility of the electrons is much larger than that of the holes, due to differences in their effective masses and their different average relaxation times [Phys. Rev. B 68, 125210 (2003)]. You can see in your calculations that convergence requires a lower momentum to converge the holes than the electrons. In these calculations, the scattering lifetimes were computed from first principles, compared to CRTA where they are constant. Phonons of momentum $\textbf{q}$ scatter the electrons, up to the defined $\textbf{q}$-point mesh. This means that increasing the $\textbf{q}$-point mesh increases the phonon modes that scatter the electrons. Using larger $\textbf{q}$-meshes will become computationally expensive, so you should aim to use the smallest $\textbf{q}$-mesh possible, while achieving convergence of $\textbf{q}$-mesh.
15.4 Questions¶
- What INCAR tag can you define the carrier densities with?
- The $\textbf{q}$-point mesh in electron-phonon calculations is the momentum of what particle?
By the end of this tutorial, you will be able to:
- Compute the band structure of Bi2Te3 with and without including spin-orbit coupling.
- Compute the different transport coefficients of Bi2Te3 in the CRTA.
- Plot the density of states, conductivity, thermopower, and ZT figure of merit as a function of chemical energy at a given temperature.
- Plot a heatmap of the ZT figure of merit as a function of carrier concentration and temperature.
16.1 Task¶
Calculate the band structure of primitive Bi2Te3 (with and without spin-polarization), then its different transport coefficients in CRTA, plotting the conductivity, DOS, Seebeck coefficient, power factor, and ZT figure of merit against the chemical potential (300 K), and produce a heatmap of the ZT figure of merit as a function of carrier concentration and temperature
A thermoelectric material is one that converts heat into electricity. Their efficiency depends on the geometry and the thermoelectric figure of merit $Z$, which is commonly expressed with temperature $T$ as $ZT$, i.e., the famous ZT figure of merit.
$$\tag{16.1} ZT=\frac{\sigma S^2 T}{\kappa_e+\kappa_l}$$
where $\sigma$ is the electrical conductivity, $S$ is the Seebeck coefficient, and $\kappa_e$ is the electronic thermal conductivity (which are ), and $\kappa_l$ is the lattice thermal conductivity, which can be calculated via the Müller-Plathe method or from phon3py.
For some context, a typical thermoelectric material has a ZT between 1 and 1.5. A new material with ZT > 2 would prove important for developing new energy harvesting techniques, accessing otherwise lost waste heat, or building compact coolers without moving parts.
The ZT formula brings together several transport coefficients: electrical conductivity, Seebeck coefficient, and thermal conductivity — each of which can be traced back to the Onsager coefficients. These, in turn, are computed from the transport distribution function, which encodes how electrons zip (or crawl) through your material.
Changing the ZT to improve thermoelectrics is a classic optimization problem [(PNAS 93, 7436 (1996)]. In this exercise, you will calculate the electronic transport coefficients of Bi2Te3 using the CRTA, in the style of Sofo and coworkers [Phys. Rev. B 68, 125210 (2003)].
16.2 Bandstructure of Bi2Te3 with and without spin-orbit coupling¶
16.2.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e16_thermoelectric_bi2te3/.
We will run the band structure calculations both without and with spin-orbit coupling (SOC) in the subdirectories bands and bands_soc, respectively. This will allow us to compare the effects of SOC on the electronic structure of Bi₂Te₃.
bands*/POSCAR
Bi2Te3
1.0
2.1932254689240596 1.2662593148768488 10.1694026073000003
-2.1932254689240596 1.2662593148768488 10.1694026073000003
-0.0000000000000001 -2.5325186297536977 10.1694026073000003
Bi Te
2 3
Direct
0.3999323159434510 0.3999323159434510 0.3999323159434510
0.6000676840565489 0.6000676840565490 0.6000676840565489
0.2098737170085619 0.2098737170085620 0.2098737170085620
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7901262829914378 0.7901262829914378 0.7901262829914378
bands/INCAR
#ground state
EDIFF = 1e-08
ISMEAR = -15
SIGMA = 0.01
LCHARG = .FALSE.
LWAVE = .FALSE.
KPAR = 4
#SOC
#LSORBIT = .TRUE.
#MAGMOM = 15*0.0
bands_soc/INCAR
#ground state
EDIFF = 1e-08
ISMEAR = -15
SIGMA = 0.01
LCHARG = .FALSE.
LWAVE = .FALSE.
KPAR = 4
#SOC
LSORBIT = .TRUE.
MAGMOM = 15*0.0
bands*/KPOINTS
Automatic generation with kspacing=0.3 div=[6, 6, 1]
0
Reduced
0.11111111111111 -0.05555555555556 -0.05555555555556
0.05555555555556 0.05555555555556 -0.11111111111111
0.33333333333333 0.33333333333333 0.33333333333333
0.00000000000000 0.00000000000000 0.00000000000000
bands*/KPOINTS_OPT
k-points for bandstructure using seekpath GAMMA-T-T-H_2-H_0-L-L-GAMMA-GAMMA-S_0-S_2-F-F-GAMMA
10
line
reciprocal
0.00000000 0.00000000 0.00000000 GAMMA
0.50000000 0.50000000 0.50000000 T
0.50000000 0.50000000 0.50000000 T
0.82299708 0.17700292 0.50000000 H_2
0.50000000 -0.17700292 0.17700292 H_0
0.50000000 0.00000000 0.00000000 L
0.50000000 0.00000000 0.00000000 L
0.00000000 0.00000000 0.00000000 GAMMA
0.00000000 0.00000000 0.00000000 GAMMA
0.33850146 -0.33850146 0.00000000 S_0
0.66149854 0.00000000 0.33850146 S_2
0.50000000 0.00000000 0.50000000 F
0.50000000 0.00000000 0.50000000 F
0.00000000 0.00000000 0.00000000 GAMMA
bands*/POTCAR
PAW Bi_d 09Feb1998
PAW Te 03Oct2001
Most of these INCAR tags you will already have come across. A few exceptions are LSORBIT and MAGMOM, which turn on spin-orbit coupling and define the initial magnetic moment for each atom, respectively. Take a look at any other tags that you do not recognise on our Wiki.
16.2.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e16_*/bands
mpirun -np 4 vasp_std
cd $TUTORIALS/electron-phonon/e16_*/bands_soc
mpirun -np 4 vasp_ncl
Once your calculations are finished, you can easily visualize and analyze the results using the py4vasp library. py4vasp offers user-friendly tools to extract, process, and plot data directly from VASP output files. With just a few lines of code, you can generate plots for band structures, density of states, and transport properties.
A typical workflow includes:
- Loading the calculation results with
py4vasp.Calculation. - Extracting key quantities such as band structure, density of states, or transport coefficients.
- Visualizing the data using built-in plotting functions, or exporting it for further analysis.
This streamlined approach simplifies post-processing and enables you to efficiently explore the electronic and transport properties of your material.
from py4vasp import Calculation
# Plot bandstructure without SOC
calc_no_soc = Calculation.from_path("e16_thermoelectric_bi2te3/bands")
band_no_soc = calc_no_soc.band.plot('kpoints_opt')
# Plot bandstructure with SOC
calc_soc = Calculation.from_path("e16_thermoelectric_bi2te3/bands_soc")
band_soc = calc_soc.band.plot('kpoints_opt')
band_no_soc.series[0].label = 'no soc'
band_soc.series[0].label = 'soc'
fig = (band_no_soc+band_soc).to_plotly()
fig.update_yaxes(range=[-1, 1])
These plots are rough due to the limited number of $\textbf{k}$-points used for the plot. Resolution could be improved by using more points. You can see that the bandgap changes when SOC is included, as well as the pocket (where the valence band rises to 0) shifting from the $\Gamma$-point without SOC to near the $T$-point with SOC, which changes the multiplicity, the effective charge, and thereby the conductivity.
16.3 Transport coefficients and ZT as a function of the chemical potential¶
16.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e16_thermoelectric_bi2te3/crta.
Bi2Te3
1.0
2.1932254689240596 1.2662593148768488 10.1694026073000003
-2.1932254689240596 1.2662593148768488 10.1694026073000003
-0.0000000000000001 -2.5325186297536977 10.1694026073000003
Bi Te
2 3
Direct
0.3999323159434510 0.3999323159434510 0.3999323159434510
0.6000676840565489 0.6000676840565490 0.6000676840565489
0.2098737170085619 0.2098737170085620 0.2098737170085620
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7901262829914378 0.7901262829914378 0.7901262829914378
#ground state
EDIFF = 1e-08
ISMEAR = -15
SIGMA = 0.01
LCHARG = .FALSE.
LREAL = .FALSE.
LSORBIT = .TRUE.
LWAVE = .FALSE.
MAGMOM = 15*0.0
#paralelization
KPAR = 4
#electron-phonon
ELPH_RUN = .TRUE.
ELPH_DRIVER = EL
ELPH_MODE=TRANSPORT
#transport
ELPH_ISMEAR = -24
EFERMI_NEDOS = 51
ELPH_SELFEN_MU_RANGE = -1.0 1.0 101
ELPH_SELFEN_TEMPS=300
ELPH_SCATTERING_APPROX=CRTA
TRANSPORT_RELAXATION_TIME = 2.2e-14
Automatic generation with kspacing=0.3 div=[6, 6, 1]
0
Reduced
0.11111111111111 -0.05555555555556 -0.05555555555556
0.05555555555556 0.05555555555556 -0.11111111111111
0.33333333333333 0.33333333333333 0.33333333333333
0.00000000000000 0.00000000000000 0.00000000000000
Automatic generation with kspacing=0.08 div=[21, 21, 3]
0
Reduced
-0.00000000000000 -0.04761904761905 -0.01587301587302
0.04761904761905 0.00000000000000 0.01587301587302
0.00000000000000 0.00000000000000 0.11111111111111
0.00000000000000 0.00000000000000 0.00000000000000
PAW Bi_d 09Feb1998
PAW Te 03Oct2001
The only new INCAR tag is ELPH_SELFEN_MU_RANGE, which sets the chemical potential as an energy shift with respect to the Fermi level, which is the limit of the chemical potential for $T \rightarrow 0$ and no additional carriers.
16.3.2 Calculation¶
You can now run your calcualtion with the following command
cd $TUTORIALS/electron-phonon/e16_*/crta
mpirun -np 4 vasp_ncl
Once it is finished you can plot the different transport coefficients using the following script:
from py4vasp import calculation, raw
import numpy as np
mus = []
conductivities = []
seebecks = []
kes = []
itemp=0
with raw.access("electron_phonon_transport",path='e16_thermoelectric_bi2te3/crta') as raw_data:
elph_transport = calculation.electron_phonon.transport.from_data(raw_data)
for iden,transport in enumerate(elph_transport):
transport_dict = transport.read_metadata()
mu = transport_dict['selfen_mu']
mus.append(mu)
transport_dict = transport.to_dict()
conductivity = transport_dict['electronic_conductivity'][itemp,:,:]
seebeck = transport_dict['seebeck'][itemp,:,:]
ke = transport_dict['electronic_thermal_conductivity'][itemp,:,:]
conductivities.append(conductivity)
seebecks.append(seebeck)
kes.append(ke)
#get density of states
import h5py
import matplotlib.pyplot as plt
vh5 = h5py.File('e16_thermoelectric_bi2te3/crta/vaspout.h5', 'r')
dos_energy = vh5['/results/electron_phonon/electrons/dos/energies'][:]
dos_total = vh5['/results/electron_phonon/electrons/dos/dos'][0,:]
efermi = vh5['/results/electron_phonon/electrons/dos/efermi'][()]
vh5.close()
mus = np.array(mus)
conductivities = np.array(conductivities)
seebecks = np.array(seebecks)
kes = np.array(kes)
#prepare plot
import matplotlib.pyplot as plt
fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, figsize=(7, 12), sharex=True)
# Conductivity subplot
ax1.set_title('Transport quantities vs chemical potential at 300K')
ax1.plot(mus, conductivities[:,0,0], label='Conductivity xx', color='#8342A4')
ax1.plot(mus, conductivities[:,2,2], '--', label='Conductivity zz', color='#35CABF')
ax1.set_ylabel('Conductivity in S/m')
ax1.legend()
# DOS subplot
ax2.plot(dos_energy-efermi,dos_total, label='DOS', color='#8342A4')
ax2.set_ylim(0,10)
ax2.set_xlim(-1,1)
ax2.set_ylabel('DOS')
ax2.legend()
# Seebeck subplot
ax3.plot(mus, seebecks[:,0,0], label='Seebeck xx', color='#8342A4')
ax3.plot(mus, seebecks[:,2,2], '--', label='Seebeck zz', color='#35CABF')
ax3.set_ylabel('Seebeck in μV/K')
ax3.legend()
# Power factor plot
ax4.plot(mus, conductivities[:,0,0]*seebecks[:,0,0]**2*1e6/100, label='Power factor xx', color='#8342A4')
ax4.plot(mus, conductivities[:,2,2]*seebecks[:,2,2]**2*1e6/100, '--', label='Power factor zz', color='#35CABF')
ax4.set_ylabel('Power factor in μW/cm K<sup>2</sup>')
ax4.legend()
# ZT plot
ax5.plot(mus, conductivities[:,0,0]*seebecks[:,0,0]**2*300/(kes[:,0,0]+1.5), label='ZT xx', color='#8342A4')
ax5.plot(mus, conductivities[:,2,2]*seebecks[:,2,2]**2*300/(kes[:,2,2]+0.7), '--', label='ZT zz', color='#35CABF')
ax5.set_xlabel('Chemical potential in eV')
ax5.set_ylabel('ZT')
ax5.legend()
plt.tight_layout()
plt.show()
All of these quantities are related by the Onsager coefficients. These figures are very similar to the ones in Fig. 6 of Sofo and coworkers [Phys. Rev. B 68, 125210 (2003)]. The agreement could be further improved by increasing the $\textbf{k}$-point mesh and using the experimental lattice parameter.
In the first plot, the conductivity $\sigma$ is effectively the DOS (second plot) smoothed out by temperature. The conductivity is greater the more carriers that there are, i.e., far from the Fermi energy.
Recollecting from the equations in the introduction of Exercise 11:
$\begin{alignedat}{2} &\text{Electrical conductivity:} &&\quad \sigma = \mathcal{L}_{11} \quad &&\text{(11.3)} \\[6pt] &\text{Seebeck coefficient:} &&\quad S = \dfrac{\mathcal{L}_{11}}{\mathcal{L}_{12}} \quad &&\text{(11.4)} \\[6pt] &\text{Electronic thermal conductivity:} &&\quad \kappa_e = \tfrac{1}{T}(\mathcal{L}_{22} - \mathcal{L}_{21}\mathcal{L}_{11}^{-1}\mathcal{L}_{12}) \quad &&\text{(11.5),} \end{alignedat}$
you can see in the third plot that the Seebeck coefficient $S$ is related to the conductivity:
$$\tag{16.2} \sigma = a(\varepsilon - \mu) + b$$
where $a$ is the slope and $b$ is the value of the conductivity at a point.
$$\tag{16.3} \sigma S^2 \propto \frac{a^2}{b^2} $$
The peaks of the Seebeck coefficient do not coincide with the peaks of the power factor. In the last plot, you can see that there are two peaks in the ZT figure of merit, one with a positive chemical potential relative to the Fermi energy and the other with a negative, i.e., the carrier peaks, holes and electrons, respectively.
In these plots, the chemical potential is shifted with respect to the Fermi level. In an experimental setting, the dopant is added to the intrinsic semiconductor to increase the number of carriers (electrons or holes).
As such, it is useful to look at the relation between chemical potential shift and carrier density. If ELPH_SELFEN_MU is specified, the number of carriers to reach that chemical shift is computed for changing temperature.
This relation is always computed: if the user specified ELPH_SELFEN_MU, then the number of carriers required to reach that the shift in chemical potential is computed, alternatively if the user specifies the number of carriers ELPH_SELFEN_CARRIER_PER_CELL or, more usually, the carrier density ELPH_SELFEN_CARRIER_DEN then the shift in chemical potential is computed:
In the plot below we read the chemical potential and the carrier density from vaspout.h5 and plot it.
from py4vasp import calculation, raw
c = calculation.from_path('e16_thermoelectric_bi2te3/crta')
chemical_potential_dict = c.electron_phonon.chemical_potential.to_dict()
efermi=chemical_potential_dict['fermi_energy']
number_of_electrons = chemical_potential_dict['carrier_density']
volume = 1e-24*c.structure.volume() #from Angstroem^3 to cm^3
carrier_density = (number_of_electrons-48) / volume
plt.plot(chemical_potential_dict['chemical_potential']-efermi, carrier_density, color='#8342A4')
plt.yscale('symlog', linthresh=1e19)
plt.ylim(-1e23, 1e23)
plt.title('Carrier density as a function of the chemical potential at 300K')
plt.ylabel(r'Carrier density in cm$^{-3}$')
plt.xlabel(r'Chemical potential $\mu$ (eV)')
plt.grid(True, alpha=0.3)
plt.show()
You can see how the carrier density and chemical potential $\mu$ are directly, though not linearly, related. All of the previous plots could therefore be plotted with respect to the carrier density, rather than the chemical potential, and they would only differ by the scaling of the x-axis. As they are not independent, you can only set the charge density (ELPH_SELFEN_CARRIER_DEN) or the chemical potential (ELPH_SELFEN_MU) in your calculation, choosing based on which property you want to measure against.
16.4 Heat map of ZT as a function of the carrier density and temperature¶
In the previous section we computed the ZT for a given temperature (300K) as a function of chemical potential shift.
Here we will plot ZT as a function of both chemical potential shift and temperature.
This is done on a heatmap with the color representing the magnitude of the figure of merit.
16.4.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron-phonon/e16_thermoelectric_bi2te3/crta_zt_map. You can generate the KPOINTS_ELPH file using the kmesh_poscar_spglib/py script, where 0.4 is the $\textbf{k}$-spacing and the name after -o is the file that it is written to:
cd $TUTORIALS/electron-phonon/e16_*/crta_zt_map/0.2_mu
python3 ../../kmesh_poscar_spglib.py -kspacing 0.4 -o KPOINTS_ELPH
0.2_mu/POSCAR and 0.2_carrier/POSCAR
Bi2Te3
1.0
2.1932254689240596 1.2662593148768488 10.1694026073000003
-2.1932254689240596 1.2662593148768488 10.1694026073000003
-0.0000000000000001 -2.5325186297536977 10.1694026073000003
Bi Te
2 3
Direct
0.3999323159434510 0.3999323159434510 0.3999323159434510
0.6000676840565489 0.6000676840565490 0.6000676840565489
0.2098737170085619 0.2098737170085620 0.2098737170085620
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7901262829914378 0.7901262829914378 0.7901262829914378
0.2_mu/INCAR
#ground state
EDIFF = 1e-08
ISMEAR = -15
SIGMA = 0.01
LCHARG = .FALSE.
LREAL = .FALSE.
LSORBIT = .TRUE.
LWAVE = .FALSE.
MAGMOM = 15*0.0
#paralelization
KPAR = 4
#electron-phonon
ELPH_RUN = .TRUE.
ELPH_DRIVER = EL
ELPH_MODE=TRANSPORT
#transport
ELPH_SELFEN_MU_RANGE = -0.5 0.5 51
ELPH_ISMEAR = -24
ELPH_SELFEN_TEMPS_RANGE=0 700 41
ELPH_SCATTERING_APPROX=CRTA
TRANSPORT_RELAXATION_TIME = 1e-14
#nedos
elph_transport_nedos = 1001
0.2_carrier/INCAR
#ground state
EDIFF = 1e-08
ISMEAR = -15
SIGMA = 0.01
LCHARG = .FALSE.
LREAL = .FALSE.
LSORBIT = .TRUE.
LWAVE = .FALSE.
MAGMOM = 15*0.0
#paralelization
KPAR = 4
#electron-phonon
ELPH_RUN = .TRUE.
ELPH_DRIVER = EL
ELPH_MODE=TRANSPORT
#transport
#ELPH_SELFEN_MU_RANGE = -1.0 1.0 51
ELPH_ISMEAR = -24
EFERMI_NEDOS = 51
ELPH_SELFEN_CARRIER_DEN_RANGE=-1e22 -1e17 51 1e17 1e22 51
ELPH_SELFEN_TEMPS_RANGE=0 700 41
ELPH_SCATTERING_APPROX=CRTA
TRANSPORT_RELAXATION_TIME = 1e-14
0.2_mu/KPOINTS and 0.2_carrier/KPOINTS
Automatic generation with kspacing=0.3 div=[6, 6, 1]
0
Reduced
0.11111111111111 -0.05555555555556 -0.05555555555556
0.05555555555556 0.05555555555556 -0.11111111111111
0.33333333333333 0.33333333333333 0.33333333333333
0.00000000000000 0.00000000000000 0.00000000000000
0.2_mu/KPOINTS_ELPH and 0.2_carrier/KPOINTS_ELPH
Automatic generation with kspacing=0.2 div=[9, 9, 2]
0
Reduced
-0.00000000000000 -0.11111111111111 -0.03703703703704
0.11111111111111 0.00000000000000 0.03703703703704
0.00000000000000 0.00000000000000 0.16666666666667
0.00000000000000 0.00000000000000 0.00000000000000
0.2_mu/POTCAR and 0.2_carrier/POTCAR
PAW Bi_d 09Feb1998
PAW Te 03Oct2001
The big difference between this INCAR file and the previous one is that instead of a single temperature ELPH_SELFEN_TEMPS=300 we specify a range of then with ELPH_SELFEN_TEMPS_RANGE=0 700 41, where the lowest temperature is given, followed by the highest, and then the number of grid points. We also specify the carrier densities that we are interested in using ELPH_SELFEN_CARRIER_DEN_RANGE, where -1e22 -1e17 51 defines that there are 51 charge densities between the upper and lower limit for holes, and 1e17 1e22 51 defines that there are 51 for electrons.
16.4.2 ZT calculation¶
You can run this calculation using:
cd $TUTORIALS/electron-phonon/e16_*/crta_zt_map/0.2_mu
mpirun -np 4 vasp_ncl
Once the calculation is finalized you can plot the results (we have set $\kappa_l$ to 1 in heat maps):
from py4vasp import calculation, raw
import numpy as np
import h5py
with raw.access("electron_phonon_transport",path='e16_thermoelectric_bi2te3/crta_zt_map/0.2_mu') as raw_data:
elph_transport = calculation.electron_phonon.transport.from_data(raw_data)
ZT = np.zeros((len(elph_transport[0].temperatures()), len(elph_transport)))
T = elph_transport[0].temperatures()
mus = []
for iden,transport in enumerate(elph_transport):
transport_dict = transport.read_metadata()
mu = transport_dict['selfen_mu']
mus.append(mu)
ZT[:,iden] = transport.figure_of_merit(kappa_lattice=1.0)
# convert dens to numpy array
muss = np.array(mus)
#get max zt
print('max ZT:',max(ZT.flatten()))
#plot zt heatmap
import matplotlib.pyplot as plt
plt.imshow(ZT, aspect='auto', cmap='viridis', extent=[min(mus), max(mus), min(T), max(T)], origin='lower')
plt.colorbar(label='ZT')
plt.xlabel('Chemical potential in eV')
plt.ylabel('Temperature in K')
plt.axhline(300, color='red', linestyle='--', label='300K')
plt.legend()
plt.show()
In the plot above, we are still computing the ZT as a function of the chemical potential, while including temperature. The previous ZT plots were performed at 300 K, where the red line is in this plot. Plotting against chemical potential is a common approach in theoretical calculations. However, in experimental settings, it is more typical to measure or control the carrier concentration directly, for example, by doping. Knowing the dopant concentration or the number of charge carriers per unit volume makes it easier to compare with experimental data and to design materials with targeted properties. Therefore, it is often more practical and insightful to plot ZT as a function of carrier density rather than chemical potential. This allows for a more direct connection between theory and experiment, and helps guide the optimization of thermoelectric materials.
We can easily perform such calculations; instead of specifying a chemical potential shift ELPH_SELFEN_MU, we can simply specify a carrier concentration ELPH_SELFEN_CARRIER_DEN. ELPH_SELFEN_CARRIER_DEN_RANGE generates a list of carrier densities for holes and electrons on a logarithmic grid between -1022 up to -1017 and then from 1017 to 1022. You could alternatively define each grid point using ELPH_SELFEN_CARRIER_DEN.
16.4.2 Carrier density calculation¶
We have already included it in the input, so you can directly run the calculation:
cd $TUTORIALS/electron-phonon/e16_*/crta_zt_map/0.2_carrier
mpirun -np 4 vasp_ncl
Once the calculation is finished, you can plot the results with the following script:
from py4vasp import calculation, raw
import numpy as np
import h5py
with raw.access("electron_phonon_transport",path='e16_thermoelectric_bi2te3/crta_zt_map/0.2_carrier') as raw_data:
elph_transport = calculation.electron_phonon.transport.from_data(raw_data)
ZT = np.zeros((len(elph_transport[0].temperatures()), len(elph_transport)))
T = elph_transport[0].temperatures()
dens = []
for iden,transport in enumerate(elph_transport):
transport_dict = transport.read_metadata()
den = transport_dict['selfen_carrier_den']
dens.append(den)
ZT[:,iden] = transport.figure_of_merit(kappa_lattice=1.0)
# convert dens to numpy array
# convert dens to log scale
dens = np.array(dens)
dens = np.sign(dens)*np.log10(np.abs(dens))
# Split dens and ZT into electron and hole parts based on sign
# Find indices for electrons (dens > 0) and holes (dens < 0)
electron_idx = dens > 0
hole_idx = dens < 0
dens_electron = dens[electron_idx]
dens_hole = dens[hole_idx]
ZT_electron = ZT[:, electron_idx]
ZT_hole = ZT[:, hole_idx]
#plot zt heatmap
import matplotlib.pyplot as plt
#get max zt
print('max hole ZT (p-doping):',max(ZT_hole.flatten()))
print('max electron ZT (n-doping):',max(ZT_electron.flatten()))
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 5), sharey=True)
# Plot for hole doping (dens < 0)
im1 = ax1.imshow(
ZT_hole,
aspect='auto',
cmap='viridis',
extent=[min(dens_hole), max(dens_hole), min(T), max(T)],
origin='lower'
)
ax1.set_title('Hole doping')
ax1.set_xlabel(r'log10(Carrier density) in cm$^{-3}$)')
ax1.set_ylabel('Temperature in K')
fig.colorbar(im1, ax=ax1, label='ZT')
# Plot for electron and hole doping
im2 = ax2.imshow(
ZT,
aspect='auto',
cmap='viridis',
extent=[np.min([dens_electron.min(), dens_hole.min()]), np.max([dens_electron.max(), dens_hole.max()]), min(T), max(T)],
origin='lower'
)
ax2.set_title('Hole + electron doping')
ax2.set_xlabel(r'log10(Carrier density) in cm$^{-3}$')
fig.colorbar(im2, ax=ax2, label='ZT')
# Plot for electron doping (dens > 0)
im3 = ax3.imshow(
ZT_electron,
aspect='auto',
cmap='viridis',
extent=[min(dens_electron), max(dens_electron), min(T), max(T)],
origin='lower'
)
ax3.set_title('Electron doping')
ax3.set_xlabel(r'log10(Carrier density) in cm$^{-3}$')
fig.colorbar(im3, ax=ax3, label='ZT')
plt.tight_layout()
plt.show()
In the earlier set of 5 plots, we visualized the Seebeck coefficient, conductivity, power factor, and ZT as a function of the chemical potential shift $\mu$ relative to the Fermi level. This is a common approach in theoretical studies, as $\mu$ is a direct output of electronic structure calculations.
In these ZT heat maps, we instead show these same transport quantities as a function of carrier density (for number of electrons or holes per unit volume), with the distribution differently spread. The key difference is the x-axis: rather than plotting against chemical potential, we plot against carrier concentration, which is the experimentally relevant parameter controlled by doping. You can see how the hole and electron carrier branches of the entire ZT heat plot are now split into two separate ones.
The relationship between chemical potential and carrier density is highly non-linear, especially near the band edges for low carrier densities. By plotting against carrier density, we provide a more intuitive and direct comparison to experimental data, and it becomes easier to identify optimal doping levels for maximizing ZT.
Despite the change in x-axis, the overall trends and features in the transport properties remain similar. Peaks and valleys in ZT, Seebeck, and conductivity as a function of $\mu$ correspond to analogous features as a function of carrier density. This demonstrates that both representations contain the same physical information, but the carrier density plot is often more practical for materials optimization and experimental comparison.
If you have time, try using a denser KPOINTS_ELPH mesh for kmesh_poscar_spglib.py in the input, e.g.:
python3 ../../kmesh_poscar_spglib.py -kspacing 0.2 -o KPOINTS_ELPH
and repeat the calculation to see how the resolution of the plots changes.
16.5 Questions¶
- What INCAR tag can you use to define a range for the chemical potential?
- How does the relationship between chemical potential and carrier density affect the interpretation of ZT plots, especially near the band edges?