by Maria Troppenz & Santiago Rigamonti for exciting neon
(Jupyter notebook by Tejus Rohatgi, Mara Voiculescu & Martin Kuban)
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.
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md
document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting
's root directory:
source tools/excitingjupyter/venv/excitingvenv/bin/activate
As a first step, you may create a running directory for the notebook.
The electronic transport coefficients, evaluated in this tutorial, are the electrical conductivity σ, 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:
\begin{equation} \tag{1} ZT=σS^2\kappa T \end{equation}
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$, i.e. $\kappa=\kappa_e+\kappa_l$. The product $σ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,
\begin{equation} \tag{2} \Xi(\epsilon) =\frac{1}{m^2_e}\sum_{n\pmb{k}} \pmb{p}_{n\pmb{k}}\pmb{p}_{n\pmb{k}}\tau_{n\pmb{k}}δ(\epsilon−\epsilon_{n\pmb{k}}) \end{equation}
It is a tensor and depends on the carrier's relaxation time $\tau_{n\pmb{k}}$the dyadic product of the charge carrier's momentum matrix elements $\pmb{p}_{n\pmb{k}}$($m_e$ 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:
\begin{equation} \tag{3} \Xi(\epsilon) =\frac{\tau}{m^2_e}\sum_{n\pmb{k}} \pmb{p}_{n\pmb{k}}\pmb{p}_{n\pmb{k}}δ(\epsilon−\epsilon_{n\pmb{k}}) \end{equation}
Using Eq. (2), the transport coefficients can be calculated as follows:
\begin{equation} \tag{4} \pmb{\sigma}(T,\mu)= e^2 \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \end{equation} \begin{equation} \tag{5} \pmb{S}(T,\mu)= \frac{ek_B}{\sigma} \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \frac{\epsilon - \mu}{k_BT} \end{equation} \begin{equation} \tag{6} \pmb{\kappa_0}(T,\mu)=k_b^2T \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \left(\frac{\epsilon - \mu}{k_BT}\right)^2 \end{equation}
Here, μ is the chemical potential, $e$ the electron charge, and $−\frac{\partial{f_0}}{\partial{\epsilon}}$ the derivative of the Fermi distribution function
\begin{equation} \tag{7} f_0(T,\mu)=\frac{1}{[exp(\epsilon−\mu)/k_BT]+1} \end{equation}
with respect to the energy $\epsilon$. $\kappa_0$ is the electronic part of the thermal conductivity for zero electrochemical potential gradient inside the sample, while the conventional electronic part of thermal conductivity, $\kappa_e$, is defined for zero electric current and 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 and its first ab initio implementation can be found in Refs. [2] and [3].
To start the tutorial, create a new, empty directory named BE-Si and move into it.
%%bash
mkdir -p BE-Si
<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 exciting_smp
with the commands
%%bash
cd BE-Si
python3 -m excitingscripts.setup.excitingroot
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
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="300">
<point coord="1.0 0.0 0.0" label="Gamma"/>
<point coord="0.625 0.375 0.0" label="K"/>
<point coord="0.5 0.5 0.0" label="X"/>
<point coord="0.0 0.0 0.0" label="Gamma"/>
<point coord="0.5 0.0 0.0" label="L"/>
</path>
</plot1d>
</bandstructure>
</properties>
and execute again.
$ time
After the calculation finished, you can visualize the band structure with the commands
$ python3 -m excitingscripts.plot.bandstructure -e -15 10
For a detailed explanation of calculating the band structure, see the tutorial Electronic band-structure and density of states.
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 Section 1) 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
do
= "fromfile"
;ngridk
= "40 40 40"
;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>
...
In order to calculate the transport coefficients, delete in your input.xml the "old" properties
block that you used for calculating the electronic band-structure:
<properties>
<bandstructure>
...
</bandstructure>
</properties>
and substitute it with
<properties>
<momentummatrix/>
<boltzequ
energyReference="efermi"
transportDfRange="-0.5 0.5"
transportDfSpacing="0.001"
temperatureRange="300 300"
temperatureSpacing="1"
chemicalPotentialRange="-0.05 0.05"
chemicalPotentialSpacing="0.0001"
evOutputEnergies="true"
siOutputUnits="true">
<etCoeffComponents>1 1</etCoeffComponents>
<etCoeffComponents>2 2</etCoeffComponents>
<etCoeffComponents>3 3</etCoeffComponents>
</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 |
---|---|---|
energyReference |
boltzequ |
If "efermi" , the energy inputs are given with respect to the Fermi energy. Else "none". |
useTransportDf |
boltzequ |
If "true" (default), the transport coefficients are calculated with the kernel function, the transport distribution function. |
transportDfRange |
boltzequ |
Frequency (or energy) range for which the transport distribution function $\Xi(\epsilon)$ (See Eq. (2) in Section 1) is calculated. |
transportDfSpacing |
boltzequ |
Spacing between the frequency (or energy) points used for the evaluation of the transport distribution function. |
temperatureRange |
boltzequ |
Temperature range in Kelvin for which the transport coefficients are calculated. |
temperatureSpacing |
boltzequ |
Spacing between temperature points in Kelvin. |
chemicalPotentialRange |
boltzequ |
Chemical potential range for which the transport coefficients are calculated. |
chemicalPotentialRange |
boltzequ |
Spacing between the chemical potential points. |
useDopingConcentration |
boltzequ |
If "true" , the chemical potential is adjusted to the doping concentration. Input parameters for the chemical potential are discarded. |
dopingConcentration |
boltzequ |
Doping concentration in cm-3. |
evOutputEnergies |
boltzequ |
If "true" , output energies are in eV. |
siOutputUnits |
boltzequ |
If "true" , output transport coefficients are given in SI units: Seebeck coefficient in V/K, electrical conductivity over relaxation time in S/(m s), thermal conductivity over relaxation time in W/(m K s). |
Using the element etCoeffComponents
, you specify the tensor component for which the transport coefficients are calculated. The element etCoeffComponents
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 exciting_smp
%%bash
cd BE-Si
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
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 range and chemical potential range defined in input.xml. They are found in the files SEEBECK_COMP.OUT
, ELECTCOND_COMP.OUT
, and THERMALCOND_COMP.OUT
, respectively, where, COMP
denotes the tensor component specified by the integer pair given in condcomp
. Furthermore, the figure of merit estimated by $Z=(\sigma S^2/\kappa_0)$ can be found in file Z_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 excitingscripts.plot.files
.
The Seebeck coefficient at the temperature of 300K can be visualized with the following commands:
%%bash
cd BE-Si
cp SEEBECK_11.OUT S11
cp SEEBECK_22.OUT S22
cp SEEBECK_33.OUT S33
python3 -m excitingscripts.plot.files -f S11 S22 S33 -cx 2 -cy 3 -ll '$S_{xx}$' '$S_{yy}$' '$S_{zz}$' -x -1 1 -lx '$\mu$ [eV]' -ly 'S [$\mu$V/K]' -ys 1e6 -rp
cd ..
The output is saved in PLOT.png
and looks as follows:
Assuming a relaxation time $\tau$ of 10 fs, the electrical conductivity can be obtained as:
%%bash
cd BE-Si
cp ELECTCOND_11.OUT E11
cp ELECTCOND_22.OUT E22
cp ELECTCOND_33.OUT E33
python3 -m excitingscripts.plot.files -f E11 E22 E33 -cx 2 -cy 3 -ll '$\sigma_{xx}$' '$\sigma_{yy}$' '$\sigma_{zz}$' -x -0.5 0.5 -y 0 9000 -lx '$\mu$ [eV]' -ly '$\sigma$ [S/cm]' -ys 1e-16 -lp 'upper center' -rp -mty 5
cd ../
Output:
And the corresponding power factor looks as:
%%bash
cd BE-Si
paste ELECTCOND_11.OUT SEEBECK_11.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11
paste ELECTCOND_22.OUT SEEBECK_22.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF22
paste ELECTCOND_33.OUT SEEBECK_33.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF33
python3 -m excitingscripts.plot.files -f PF11 PF22 PF33 -cx 2 -cy 3 -ll '${\sigma S^2}_{xx}$' '${\sigma S^2}_{yy}$' '${\sigma S^2}_{zz}$' -x -0.5 0.5 -y 0 20 -lx '$\mu$ [eV]' -ly '$\sigma S^2$ [$\mu W\,cm^{-1}\,K^{-2}$]' -ys 1e-10 -lp 'upper center' -rp
cd ../
Output:
In the following, we look at the temperature dependence of transport coefficients 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
do
to "skip" in the element groundstate
The element properties
of your input.xml can be changed to
<properties>
<!--momentummatrix/-->
<boltzequ
energyReference="efermi"
transportDfRange="-0.5 0.5"
transportDfSpacing="0.001"
temperatureRange="300 1200"
temperatureSpacing="10"
chemicalPotentialRange="-0.005 0.005"
chemicalPotentialSpacing="0.005"
evOutputEnergies="true"
siOutputUnits="true">
<etCoeffComponents>1 1</etCoeffComponents>
</boltzequ>
</properties>
<properties>
<momentummatrix/>
<boltzequ
energyReference="efermi"
transportDfRange="-0.5 0.5"
transportDfSpacing="0.001"
temperatureRange="300 300"
temperatureSpacing="1"
chemicalPotentialRange="-0.05 0.05"
chemicalPotentialSpacing="0.0001"
evOutputEnergies="true"
siOutputUnits="true">
<etCoeffComponents>1 1</etCoeffComponents>
<etCoeffComponents>2 2</etCoeffComponents>
<etCoeffComponents>3 3</etCoeffComponents>
</boltzequ>
</properties>
The transport coefficients are calculated from 300 to 1200 K and for the chemical potentials $\mu = -0.136$ eV, 0 eV, and $0.136$ eV with respect to the Fermi energy. Here, we compute only the 1 1
-component of the transport coefficient's tensor. With the modified input.xml, we can rerun exciting_smp
%%bash
cd BE-Si
time $EXCITINGROOT/bin/exciting_smp &
cd ../
The Seebeck coefficient for $\mu = -0.136$ eV and $\mu = 0.136$ eV can be visualized as:
%%bash
cd BE-Si
grep ' -0.136' SEEBECK_11.OUT > S11M
grep ' 0.136' SEEBECK_11.OUT > S11P
python3 -m excitingscripts.plot.files -f S11P S11M -cy 3 -ll '$S_{xx}$@$\mu$= 0.136eV' '$S_{xx}$@$\mu$=-0.136eV' -x 300 1200 -lx 'T [K]' -ly 'S [$\mu$V/K]' -ys 1e6 -lp 5 -rp
cd ../
Output
Again, using τ = 10 fs, the electrical conductivity for these μ's is obtained as:
%%bash
cd BE-Si
grep ' -0.136' ELECTCOND_11.OUT > E11M
grep ' 0.136' ELECTCOND_11.OUT > E11P
python3 -m excitingscripts.plot.files -f E11P E11M -cy 3 -ll '$\sigma_{xx}$@$\mu$= 0.136eV' '$\sigma_{xx}$@$\mu$=-0.136eV' -x 300 1200 -lx 'T [K]' -ly '$\sigma$ [S/cm]' -ys 1e-16 -rp
cd ../
Output:
And the corresponding power factor looks as:
%%bash
cd BE-Si
paste ELECTCOND_11.OUT SEEBECK_11.OUT | grep ' -0.136' | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11M
paste ELECTCOND_11.OUT SEEBECK_11.OUT | grep ' 0.136' | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11P
python3 -m excitingscripts.plot.files -f PF11P PF11M -cy 3 -ll '${\sigma S^2}_{xx}$@$\mu$= 0.136eV' '${\sigma S^2}_{xx}$@$\mu$=-0.136eV' -x 300 1200 -lx 'T [K]' -ly '$\sigma S^2$ [$\mu W\,cm^{-1}\,K^{-2}$]' -ys 1e-10 -rp
cd ../
Output: