by Maria Troppenz, Patrick Dieu, & Santiago Rigamonti for exciting nitrogen
Purpose: In this tutorial you will learn how to calculate electronic transport coefficients with exciting. The coefficients are obtained by solving the linearized Boltzmann equation in the constant relaxation time approximation. They are calculated for bulk silicon.
0. Theoretical background
The electronic transport coefficients, evaluated in this tutorial, are the electrical conductivity $\sigma$, the Seebeck coefficient $S$, and the electronic part of the thermal conductivity $\kappa_e$. They are related to the electronic structure of the material and directly affect the material's thermoelectric efficiency through the dimensionless figure of merit:
(1)
\begin{align} Z T= \frac { \sigma S^2}{\kappa} T ~. \end{align}
Here, $T$ is the temperature and $\kappa$ the total thermal conductivity, which is given by the sum of electronic contribution $\kappa_e$ and the lattice contribution $\kappa_l$. The product $\sigma S^2$ is also called the power factor.
The transport coefficients are obtained by solving the linearized Boltzmann's equation in the relaxation time approximation (for further information see Ref.1). The kernel of the transport coefficients is the transport distribution function,
(2)
\begin{align} \pmb{\Xi} (\epsilon)= \frac{1}{m_e^2} \sum_{ n \pmb{k} }\, { \pmb{p}}_{n \pmb{k}} \, { \pmb{p} }_{n \pmb{k}} \, \tau_{ n \pmb{k} } \, \delta(\epsilon-\epsilon_{ n \pmb{k} })~. \end{align}
It is a tensor and depends on the carrier's relaxation time $\tau_{n \pmb{k}}$ and the dyadic product of the charge carrier's momentum matrix elements ${ \pmb{p} }_{n \pmb{k} }$ ($m_e$ is the electron mass). The sum is over all states with the band index $n$ and the wave vector $\pmb{k}$ in the first Brillouin zone. In the constant relaxation time approximation, $\tau$ is independent from $n$ and $\pmb{k}$ and the transport distribution function is $\pmb{\Xi} (\epsilon)=\tau \frac{1}{m_e^2} \sum_{ n \pmb{k} }\, { \pmb{p}}_{n \pmb{k}} \, { \pmb{p} }_{n \pmb{k}} \, \delta(\epsilon-\epsilon_{ n \pmb{k} })$.
Using Eq.(2), the transport coefficients can be calculated as follows:
(3)
\begin{align} \pmb{\sigma} (T,\mu) = e^2 \int\! \mathrm{d} \epsilon \, \left( - \frac{\partial f_{0}(T,\mu)} { \partial \epsilon } \right) \, \pmb{\Xi}(\epsilon) \end{align}
(4)
\begin{align} \pmb{S} (T,\mu) = \frac{ e\, k_B}{\sigma} \int\! \mathrm{d} \epsilon \, \left( - \frac{\partial f_{0}(T,\mu)} { \partial \epsilon } \right) \, \pmb{\Xi}(\epsilon) \, \frac{\epsilon - \mu}{k_B T} \end{align}
(5)
\begin{align} \pmb{\kappa}_0 (T, \mu) = k_B^2 T \int\! \mathrm{d} \epsilon \, \left( - \frac{\partial f_{0}(T,\mu)} { \partial \epsilon } \right) \, \pmb{\Xi}(\epsilon) \, \left(\frac{\epsilon - \mu}{k_B T} \right)^2 \end{align}
Here, $\mu$ is the chemical potential, $e$ the electron charge, and $- \partial f_0 / \partial \epsilon$ the derivative of the Fermi distribution function $f_0 (T, \mu) = [ (\exp( (\epsilon-\mu) / k_B T) +1 ]^{-1}$ with respect to the energy $\epsilon$. Using the equations above, the electronic part of thermal conductivity for zero electric current $\kappa_e$ is obtained by the relation $\kappa_e = \kappa_0 - T \sigma S^2$. The lattice part of the thermal conductivity still remains to be calculated by other methods.
Details about the method can be found in Ref.2 and 3.
1. Define relevant environment variables
Read the following paragraphs before starting with the rest of this tutorial!
Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. The following script is used in this tutorial:
- PLOT-transport.py: Python script for visualizating the transport coefficients.
From now on the symbol $ will indicate the shell prompt.
Important note: All input parameters are in atomic units!
2. Transport coefficients of Bulk Silicon
To start the tutorial, create a new, empty directory named BE-Si and move into it.
i) Preparation of the calculation
Save the following lines as input.xml:
<input>
<title>Bulk Silicon</title>
<structure speciespath="$EXCITINGROOT/species/">
<crystal scale="10.26">
<basevect>0.0 0.5 0.5</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.5 0.5 0.0</basevect>
</crystal>
<species speciesfile="Si.xml">
<atom coord="0.00 0.00 0.00"></atom>
<atom coord="0.25 0.25 0.25"></atom>
</species>
</structure>
<groundstate
do="fromscratch"
xctype="GGA_PBE_SOL"
ngridk="4 4 4"
rgkmax="8.0"
gmaxvr="14.0"
maxscl="100">
</groundstate>
</input>
Set up $EXCITINGROOT to the correct root directory in speciespath and execute excitingser with the commands
$ SETUP-excitingroot.sh
$ time excitingser
To gain a better understanding of the later calculated transport coefficients of silicon, you may want to look at the band structure. For this, please expand the link below and follow the instructions.
Add the element properties in your input.xml with the following lines
...
<properties>
<bandstructure>
<plot1d>
<path steps="100">
<point coord="0.250 0.750 0.500" label="W" />
<point coord="0.000 0.000 0.000" label="GAMMA"/>
<point coord="0.000 0.500 0.500" label="X" />
<point coord="0.250 0.750 0.500" label="W"/>
<point coord="0.000 0.500 0.000" label="L" />
<point coord="0.000 0.000 0.000" label="GAMMA"/>
</path>
</plot1d>
</bandstructure>
</properties>
...
and execute
excitingser again.
After the calculation finished, you can visualize the band structure with the commands
$ xsltproc $EXCITINGVISUAL/xmlband2agr.xsl bandstructure.xml
$ xmgrace Bulk\ Silicon_bandstructure.agr >/dev/null 2>&1
For a detailed explanation of calculating the band structure, see tutorial nitrogen-electronic-structure-calculations.
ii) Calculation of the transport coefficients at room temperature
In the following, we calculate the transport coefficients with respect to the chemical potential at a constant temperature (here 300K). The equations of the transport coefficients (see 0.Theoretical Background) contain the derivative of the Fermi distribution function $f_0$. That means that only k-points with eigenenergies near the Fermi surface contribute to the coefficients. In order to achieve a sufficient representation of the energy bands in this area, a fine k-grid is needed.
To obtain eigenenergies for a fine k-grid, we perform a single self-consistency loop on top of the previous ground-state calculation by changing in input.xml
- the attribute do = "fromfile";
- the attribute ngridk = "40 40 40";
- the attribute maxscl = "1";
in the element groundstate (for more details see Input Reference).
The element groundstate in input.xml should now look like this:
...
<groundstate
do="fromfile"
xctype="GGA_PBE_SOL"
ngridk="40 40 40"
rgkmax="8.0"
gmaxvr="14.0"
maxscl="1">
</groundstate>
...
The transport coefficients are calculated by adding the following lines in the element properties of your input.xml
...
<properties>
<momentummatrix/>
<boltzequ
windtdf="-1.0 0.5"
nwtdf="1500"
windtemp="300 300"
tgrid="1"
windmu="0.1691 0.2491"
mugrid="200"
tevout="true"
tsiout="true">
<condcomp>1 1</condcomp>
<condcomp>2 2</condcomp>
<condcomp>3 3</condcomp>
</boltzequ>
</properties>
...
With inserting the element momentummatrix, the momentum matrix elements are calculated. Inside the element boltzequ, you define the attributes needed for the calculation of the transport coefficients.
An explanation of the attributes can be found by expanding the link below or in Input Reference.
| Attribute |
Element |
Description |
| windtdf |
boltzequ |
Energy window in which the transport distribution function $\pmb{\Xi}(\epsilon)$ (Eq.(2)) is calculated. |
| nwtdf |
boltzequ |
Number of energy points in the energy window for the transport distribution function. |
| windtemp |
boltzequ |
Temperature range in which the transport coefficients are calculated. |
| tgrid |
boltzequ |
Number of temperature values in the temperature range. |
| windmu |
boltzequ |
Range of the chemical potential in which the transport coefficients are calculated. |
| mugrid |
boltzequ |
Number of values in the range of the chemical potential. |
| tevout |
boltzequ |
If "true", energy outputs in the output files for transport are in eV. |
| tsiout |
boltzequ |
If "true", the transport coefficients in the output files are given in SI units: Seebeck coefficient in V/K , electrical conductivity over tau in S/(m s), thermal conductivity over tau in W/(mK s). |
|
Using the element
condcomp, you specify the tensor component for which the transport coefficients are calculated. The element
condcomp needs, as input, an integer pair,
e.g.,
1 1 stands for the first diagonal component of the tensor (
1 =
x-direction,
2 =
y-direction, and
3 =
z-direction).
By executing excitingser
the Seebeck coefficient
$S$, the electrical conductivity over the relaxation time
$\sigma /\tau$, and the thermal conductivity over the relaxation time
$\kappa_0 / \tau$ are calculated for the temperature values and chemical potential values defined in
input.xml. The calculated transport coefficients are in the files
SEEBECK_COMP.OUT,
ELECTCOND_COMP.OUT, and
THERMALCOND_COMP.OUT, respectively, while
COMP denotes the tensor component specified by the integer pair given in
condcomp. Furthermore, the figure of merit estimated by
$ZT = ( \sigma S^2/\kappa_0) \,T$ can be found in file
ZT_COMP.OUT.
Please note: $\sigma$ and $\kappa_0$ are not directly given, but the output files contain the coefficients $\sigma /\tau$ and $\kappa_0 / \tau$. $\tau$ itself is constant (does not depend on electronic structure with $\pmb{k}$ and $n$), and can be multiplied afterwards.
Now, the transport coefficients are calculated and you can plot them by using PLOT-transport.py.
The Seebeck coefficient, the electrical conductivity over $\tau$ and the power factor for the temperature of 300K can be visualized with the following commands:
$ PLOT-transport.py single[S] T=300.0K
$ PLOT-transport.py single[sigma] T=300.0K
$ PLOT-transport.py single[PF] T=300.0K,t=1e-14s
The output looks as follows:
A detailed description of the plot script is given in the link below.
ii.i) Description how to visualize transport coefficients
To visualize one of the transport coefficients at a constant temperature, run the command:
$ PLOT-transport single[<X>] T=<Y>K
At position
<X>, you put either the symbol
S for the Seebeck coefficient,
sigma for the electrical conductivity,
PF for the power factor
$\sigma S^2$,
kappa for the thermal conductivity, or
ZT for the figure of merit. At position
<Y>, the temperature is defined. The value of the temperature, you get from the first column in output files of the transport properties
To visualize one of the transport coefficients at a constant chemical potential, run
$ PLOT-transport single[<X>] mu=<Z>eV
or
$ PLOT-transport single[<X>] mu=<Z>Ha
depending on the attribute
tevout in your input file (eV if
tevout =
"true" and Ha if
tevout =
"false"). At position
<Z>, now, the value of the chemical potential is defined and can be extracted from the second column of the output files of the transport properties.
So far, the plots show only the coefficients over the relaxation time $\tau$. You can define a relaxation time <TAU>, which usually is in the order of 10 fs $\approx 1 \cdot 10^{-14}$, with using one of these commands:
$ PLOT-transport single[<X>] T=<Y>K,t=<TAU>s
$ PLOT-transport single[<X>] mu=<Z>eV,t=<TAU>s
$ PLOT-transport single[<X>] mu=<Z>Ha,t=<TAU>s
Advanced users can plot many transport coefficients at the same time by, e.g., using the command:
$ PLOT-transport pattern[S,sigma,PF,ZT] T=<Y>K
The brackets after "pattern" can contain any of the transport coefficients, separated each by a comma.
iii) Calculation of the transport coefficients at different doping concentrations
In the following, we look at the transport coefficients with respect to the temperature at a constant chemical potential.
We can use the fine k-grid obtained in part ii) and skip the ground-state calculation by changing in input.xml
- the attribute do = "skip";
in the element groundstate. Moreover, the momentum matrix elements does not need to be calculated again and, e.g., can be removed or commented.
The element properties of your input.xml can be changed to
...
<properties>
<boltzequ
windtdf="-1.0 0.5"
nwtdf="1500"
windtemp="300 1200"
tgrid="90"
windmu="0.2026 0.2107"
mugrid="3"
tevout="true"
tsiout="true">
<condcomp>1 1</condcomp>
<condcomp>2 2</condcomp>
<condcomp>3 3</condcomp>
</boltzequ>
</properties>
...
With the modified input.xml, we can rerun excitingser
The Fermi level of the ground state for silicon is at around 5.6231 eV and represents the charge neutrality level (value can be extracted from EFERMI.OUT). The chemical potential can be set to values below this energy (hole-doped silicon), or to values above this energy (electron-doped silicon), in order to simulate different doping concentrations for bulk silicon.
In the transport output files, there are different values of the chemical potential.
We can look at some transport coefficients for the hole-doped silicon with these commands:
$ PLOT-transport.py single[S] mu=5.6231eV
$ PLOT-transport.py single[sigma] mu=5.6231eV
$ PLOT-transport.py single[PF] mu=5.6231eV,t=1e-14s
The output looks as follows:
Results for transport coefficients of silicon alloys calculated with Boltzmann's equation are published in Ref.4, and can be compare to your calculation.
Bibliography
1. Method described in detail in: B. R. Nag, "Electron Transport in Compound Semiconductors" (Springer Verlag, Berlin, 1996).
2. G. D. Mahan and J. O. Sofo, "The best thermoelectric", Proc. Natl. Aad. Sci. USA 93, 7436 (1996).
3. T. J. Scheidemantel, C. Ambrosch-Draxl, T. Thonhauser, H. V. Badding, and J. O. Sofo, "Transport coefficients from first-principles calculations", Phys. Rev. B 68, 125210 (2003).
4. S. Nakajima, J. Phys. Chem. Solids 24, 479 (1963).
