Brillouin Zone Integration: Linear Tetrahedron Method
Post
Cancel

# Brillouin Zone Integration: Linear Tetrahedron Method

## Introduction

The translational symmetry of solids resulits in a quantum number, i.e. the crystal momentum $\mathbf{k}$. Consequently many quantities of the crystal, e.g. total energy and density of states, require integration over the Brillouin Zone (BZ). The integrations are generally of the form

$$$\label{eq:bz_int} \frac{1}{\Omega_{BZ}}\int_{BZ} \,M_{\mathbf{k}}\cdot f(\varepsilon_\mathbf{k})\,\mathrm{d}\mathbf{k}$$$

where $M_{\mathbf{k}}$ is the matrix element and $f(\varepsilon_\mathbf{k})$ is Heaviside step function $\Theta(E - \varepsilon_\mathbf{k})$ or Dirac function $\delta(E - \varepsilon_\mathbf{k})$.

There are generally two categories of methods to perform the integration: the special-point method and the tetrahedron method, where the former method is often combined with smearing method to improve convergence.1 The most widely used sets of special points are those of Monkhorst and Pack,2 however, it is beyond the scope of the current post to discuss it.

## Tetrahedron Method

• In the tetrahedron method, the reciprocal space is first divided into small sub-meshes, which are then further divided into tetrahedra, as is shown in Figure 1.

Figure 1. Breakup of a sub-mesh into six tetrahedra. Generated by Plotly, another script using matplotlib can be found here.
• Within each tetrahedron, the quantity $M_{\mathbf{k}}$ and $\varepsilon_\mathbf{k}$ are linearly or quadraticly interpolated and as a result, the BZ integration can be done analytically. We will focus on the linear interpolation case, known as linear tetrahedron method.3

### Linear Tetrahedron method

In the linear tetrahedron method, the quantities, $M_\mathrm{k}$ and $\varepsilon_\mathbf{k}$, are linearly interpolated within the tetrahedron, i.e.

$$$\varepsilon(x, y, z) = a + b\cdot x + c\cdot y + d\cdot z$$$

where $x, y, z$ are the three component of $\mathbf{k}$. With the values at the four vertices, e.g. $\varepsilon_i (i=1,4)$, the coeffiecients $a$, $b$, $c$ and $d$ can be readily obtained. One can also use the Barycentric coordinates4 of the tetrahedron and express the quantity as

$$$\label{eq:barycentric_exp} \varepsilon(e, u, v) = \varepsilon_1\cdot(1 - e - u - v) + \varepsilon_2 \cdot e + \varepsilon_3 \cdot u + \varepsilon_4 \cdot v$$$

where

\begin{align} x &= x_1\cdot(1 - e - u - v) + x_2 \cdot e + x_3 \cdot u + x_4 \cdot v \\ y &= y_1\cdot(1 - e - u - v) + y_2 \cdot e + y_3 \cdot u + y_4 \cdot v \\ z &= z_1\cdot(1 - e - u - v) + z_2 \cdot e + z_3 \cdot u + z_4 \cdot v \\ \end{align}

and $e, u, v \in [0, 1]$.

Now the contribution within the tetrahedron to the BZ integration becomes

$$$\frac{1}{\Omega_{BZ}}\int_{BZ}\, M(e, u, v) \cdot f(\varepsilon(e, u, v)) \cdot \frac{\partial(x, y, z)}{\partial(e, u, v)} \, \mathrm{d}e \mathrm{d}u \mathrm{d}v$$$

where $\frac{\partial(x, y, z)}{\partial(e, u, v)}$ is the Jacobi determinant and one can readily show that it equals to $6\Omega_T$ where $\Omega_T$ is the volume of the tetrahedron.5

### Density of States (DOS)

For the density of state (DOS), the quantity $M_\mathbf{k} =1$ and $f(\varepsilon_\mathbf{k})$ is the Dirac function. In this case, the contribution of the i-th tetrahedron $T_i$ to the DOS is

\begin{align} \label{eq:dos} \rho(E) &= \frac{1}{\Omega_{BZ}}\int_{T_i}\, \delta(E - \varepsilon(e, u, v)) \cdot \frac{\partial(x, y, z)}{\partial(e, u, v)} \, \mathrm{d}e \mathrm{d}u \mathrm{d}v \\[9pt] &= \frac{6\Omega_T}{\Omega_{BZ}}\int_{T_i}\, \frac{1}{|\nabla \varepsilon(e, u, v)|} \, \mathrm{d}S\Bigr|_{\varepsilon = E} \end{align}

where the volume integration over the BZ becomes a surface integration on an iso-value plane. Moreover, let us assumed $\varepsilon_i$ is ordered according to increasing values, i.e. $\varepsilon_1 < \varepsilon_2 < \varepsilon_3 < \varepsilon_4$ and denote $\varepsilon_{ji} = \varepsilon_j - \varepsilon_i$, where $\varepsilon_0 = E$.

The norm of the derivative can be reaily get from Eq.\eqref{eq:barycentric_exp}

$$$|\nabla\varepsilon(e, u, v)| = \sqrt{\varepsilon_{21}^2+\varepsilon_{31}^2+\varepsilon_{41}^2}$$$

Figure 2. Isosurfaces of different energy $E$ shown in red plane. Generated by Plotly.

For different value of the energy $E$, the iso-value surfaces intersected with the tetrahedron are differnt:

• For $\varepsilon_4 < E$ or $E > \varepsilon_E$

No intersection between the iso-surface and the tetrahedron, hence no contribution to DOS within this tetrahedron.

• For $\varepsilon_1 < E < \varepsilon_2$

The iso-surface is a triangle, as is shown in Figure 2. The coordinates of vertices of the triangle $\vec{c_i}$ can be readily obtained from Eq.\eqref{eq:barycentric_exp}

\begin{align} \vec{c_1} \quad&\Rightarrow\quad \begin{pmatrix} \dfrac{\varepsilon_{01}}{\varepsilon_{21}}, & 0, & 0 \end{pmatrix} \\[9pt] \vec{c_2} \quad&\Rightarrow\quad \begin{pmatrix} 0, & \dfrac{\varepsilon_{01}}{\varepsilon_{31}}, & 0 \end{pmatrix} \\[9pt] \vec{c_3} \quad&\Rightarrow\quad \begin{pmatrix} 0, & 0, & \dfrac{\varepsilon_{01}}{\varepsilon_{41}} \end{pmatrix} \\[9pt] \end{align}

The area of the iso-surface is then given by

$$$\frac{1}{2}\cdot|(\vec{c_1} - \vec{c_3})\times(\vec{c_2} - \vec{c_3})| = \frac{1}{2} \cdot \varepsilon_{01}^2 \cdot \frac{\sqrt{\varepsilon_{21}^2+\varepsilon_{31}^2+\varepsilon_{41}^2}} {\varepsilon_{21}\varepsilon_{31}\varepsilon_{41}}$$$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 #!/usr/bin/env python # -*- coding: utf-8 -*- import sympy as sp e01, e21, e31, e41 = sp.symbols('e01, e21, e31, e41') c1 = sp.Matrix([e01/e21, 0, 0]) c2 = sp.Matrix([0, e01/e31, 0]) c3 = sp.Matrix([0, 0, e01/e41]) s1 = (c1 - c3).cross(c2 - c3) sp.pprint(sp.simplify(sp.sqrt(s1.dot(s1)))) # ___________________________ # ╱ 4 ⎛ 2 2 2⎞ # ╱ e₀₁ ⋅⎝e₂₁ + e₃₁ + e₄₁ ⎠ # ╱ ───────────────────────── # ╱ 2 2 2 # ╲╱ e₂₁ ⋅e₃₁ ⋅e₄₁ # 

Therefore, the contribution to the DOS is

$$$\rho_1(E) = \frac{\Omega_T}{\Omega_{BZ}} \cdot \frac{3\cdot \varepsilon_{01}^2} {\varepsilon_{21}\varepsilon_{31}\varepsilon_{41}}$$$
• For $\varepsilon_2 < E < \varepsilon_3$

The iso-surface is a tetragon, as is shown in Figure 2. Again, the coordinates of vertices of the tetragon are given by

\begin{align} \vec{c_1} \quad&\Rightarrow\quad \begin{pmatrix} \dfrac{\varepsilon_{40}}{\varepsilon_{42}}, & 0, & -\dfrac{\varepsilon_{20}}{\varepsilon_{42}} \end{pmatrix} \\[9pt] \vec{c_2} \quad&\Rightarrow\quad \begin{pmatrix} \dfrac{\varepsilon_{30}}{\varepsilon_{32}},& -\dfrac{\varepsilon_{20}}{\varepsilon_{32}},& 0 \end{pmatrix} \\[9pt] \vec{c_3} \quad&\Rightarrow\quad \begin{pmatrix} 0,& \dfrac{\varepsilon_{01}}{\varepsilon_{31}},& 0 \end{pmatrix} \\[9pt] \vec{c_4} \quad&\Rightarrow\quad \begin{pmatrix} 0,& 0,& \dfrac{\varepsilon_{01}}{\varepsilon_{41}} \end{pmatrix} \\[9pt] \end{align}

The area of the iso-surface is then given by

\begin{align*} S &= \frac{1}{2}\cdot \Bigl| \left[ (\vec{c_1} - \vec{c_4}) + (\vec{c_2} - \vec{c_4}) \right] \times(\vec{c_2} - \vec{c_4}) \Bigr| \\[9pt] &= \frac{1}{2} \cdot \frac{\sqrt{\varepsilon_{21}^2+\varepsilon_{31}^2+\varepsilon_{41}^2}} {\varepsilon_{31}\varepsilon_{41}} \cdot \left[ \varepsilon_{21} - 2\varepsilon_{20} - \frac{\varepsilon_{20}^2\cdot(\varepsilon_{31} + \varepsilon_{42})} {\varepsilon_{32}\varepsilon_{42}} \right] \end{align*}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 #!/usr/bin/env python # -*- coding: utf-8 -*- import sympy as sp e43, e42, e41, e40, e32, e31, e30, e21, e20, e10 = sp.symbols( 'e43, e42, e41, e40, e32, e31, e30, e21, e20, e10' ) c1 = sp.Matrix([e40 / e42, 0, - e20 / e42]) c2 = sp.Matrix([e30 / e32, -e20 / e32, 0]) c3 = sp.Matrix([0, -e10 / e31, 0]) c4 = sp.Matrix([0, 0, -e10 / e41]) s1 = (c1 + c3 - 2 * c4).cross(c2 - c4) sp.pprint( sp.simplify( sp.sqrt(s1.dot(s1)) ) ) # _______________________________________________________________________________________________________________________________________ # ╱ 2 # ╱ 2 2 2 2 ⎛ 2 ⎞ # ╱ e₃₁ ⋅(e₁₀⋅e₃₂⋅e₄₀ - e₃₀⋅(2⋅e₁₀⋅e₄₂ - e₂₀⋅e₄₁)) + e₄₁ ⋅(e₁₀⋅e₃₀⋅e₄₂ - e₂₀⋅e₃₁⋅e₄₀) + ⎝e₁₀ ⋅e₃₂⋅e₄₂ - e₂₀⋅e₃₁⋅(2⋅e₁₀⋅e₄₂ - e₂₀⋅e₄₁)⎠ # ╱ ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── # ╱ 2 2 2 2 # ╲╱ e₃₁ ⋅e₃₂ ⋅e₄₁ ⋅e₄₂ # one can show that the expression within the braces are identical, i.e. # 2 2 # (e₁₀⋅e₃₂⋅e₄₀ - e₃₀⋅(2⋅e₁₀⋅e₄₂ - e₂₀⋅e₄₁)) = (e₁₀⋅e₃₀⋅e₄₂ - e₂₀⋅e₃₁⋅e₄₀) # 2 # ⎛ 2 ⎞ 2 2 # ⎝e₁₀ ⋅e₃₂⋅e₄₂ - e₂₀⋅e₃₁⋅(2⋅e₁₀⋅e₄₂ - e₂₀⋅e₄₁)⎠ = e₂₁ ⋅(e₀₁⋅e₃₀⋅e₄₂ + e₂₀⋅e₃₁⋅e₄₀) # sp.pprint( sp.expand( (e42*(e32 + e20)*(e21 - e20) - e20*e31*(e42 + e20)) ) ) # 0 < (e₀₁⋅e₃₀⋅e₄₂ + e₂₀⋅e₃₁⋅e₄₀) = # # 2 2 # - e₂₀ ⋅e₃₁ - e₂₀ ⋅e₄₂ + e₂₀⋅e₂₁⋅e₄₂ - e₂₀⋅e₃₁⋅e₄₂ - e₂₀⋅e₃₂⋅e₄₂ + e₂₁⋅e₃₂⋅e₄₂ = # # 2 # - e₂₀ ⋅(e₃₁ + e₄₂) - 2⋅e₂₀⋅e₃₂⋅e₄₂ + e₂₁⋅e₃₂⋅e₄₂ 

Therefore, the contribution to the DOS is

$$$\rho_2(E) = \frac{\Omega_T}{\Omega_{BZ}} \cdot \frac{1}{\varepsilon_{31}\varepsilon_{41}} \cdot \left[ 3\varepsilon_{21} - 6\varepsilon_{20} - 3\frac{\varepsilon_{20}^2\cdot(\varepsilon_{31} + \varepsilon_{42})} {\varepsilon_{32}\varepsilon_{42}} \right]$$$
• For $\varepsilon_3 < E < \varepsilon_4$

The iso-surface is a triangle, as is shown in Figure 2. The coordinates of vertices of the triangle $\vec{c_i}$

\begin{align} \vec{c_1} \quad&\Rightarrow\quad \begin{pmatrix} \dfrac{\varepsilon_{40}}{\varepsilon_{42}}, & 0, & -\dfrac{\varepsilon_{20}}{\varepsilon_{42}} \end{pmatrix} \\[9pt] \vec{c_2} \quad&\Rightarrow\quad \begin{pmatrix} 0, & \dfrac{\varepsilon_{40}}{\varepsilon_{43}}, & -\dfrac{\varepsilon_{30}}{\varepsilon_{43}} \end{pmatrix} \\[9pt] \vec{c_3} \quad&\Rightarrow\quad \begin{pmatrix} 0, & 0, & \dfrac{\varepsilon_{01}}{\varepsilon_{41}} \end{pmatrix} \\[9pt] \end{align}

The area of the iso-surface is then given by

$$$\frac{1}{2}\cdot|(\vec{c_1} - \vec{c_3})\times(\vec{c_2} - \vec{c_3})| = \frac{1}{2} \cdot \varepsilon_{40}^2 \cdot \frac{\sqrt{\varepsilon_{21}^2+\varepsilon_{31}^2+\varepsilon_{41}^2}} {\varepsilon_{41}\varepsilon_{42}\varepsilon_{43}}$$$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 #!/usr/bin/env python # -*- coding: utf-8 -*- import sympy as sp e43, e42, e41, e40, e30, e20, e10 = sp.symbols( 'e43, e42, e41, e40, e30, e20, e10' ) c1 = sp.Matrix([e40 / e42, 0, - e20 / e42]) c2 = sp.Matrix([0, e40 / e43, - e30 / e43]) c3 = sp.Matrix([0, 0, -e10 / e41]) s1 = (c1 - c3).cross(c2 - c3) sp.pprint(sp.simplify(sp.sqrt(s1.dot(s1)))) # ________________________________________________________________ # ╱ 2 ⎛ 2 2 2 2⎞ # ╱ e₄₀ ⋅⎝e₄₀ ⋅e₄₁ + (e₁₀⋅e₄₂ - e₂₀⋅e₄₁) + (e₁₀⋅e₄₃ - e₃₀⋅e₄₁) ⎠ # ╱ ────────────────────────────────────────────────────────────── # ╱ 2 2 2 # ╲╱ e₄₁ ⋅e₄₂ ⋅e₄₃ # one can easily show that # (e₁₀⋅e₄₂ - e₂₀⋅e₄₁) = -e₄₀⋅e₂₁ # and # (e₁₀⋅e₄₃ - e₃₀⋅e₄₁) = -e₄₀⋅e₃₁ 

Therefore, the contribution to the DOS is

$$$\rho_3(E) = \frac{\Omega_T}{\Omega_{BZ}} \cdot \frac{3\cdot\varepsilon_{40}^2} {\varepsilon_{41}\varepsilon_{42}\varepsilon_{43}}$$$

### No. of States

The number of states $n(E)$ is simply the integration of the DOS, i.e.

$$$n(E) = \int_0^E \rho(E)\, \mathrm{d}E$$$

As a result, the contributions of the i-th tetrahedron to the no. of states are

• For $E < \varepsilon_1$
$$$n_0(E) = 0$$$
• For $\varepsilon_1 < E \le \varepsilon_2$
$$$n_1(E) = \int_{\varepsilon_1}^E \, \rho_1(E)\, \mathrm{d}E = \frac{\Omega_T}{\Omega_{BZ}} \cdot \frac{\varepsilon_{01}^3} {\varepsilon_{21}\varepsilon_{31}\varepsilon_{41}}$$$
• For $\varepsilon_2 < E \le \varepsilon_3$
\begin{align} n_2(E) &= n_1(\varepsilon_2) + \int_{\varepsilon_2}^E \, \rho_2(E)\, \mathrm{d}E \nonumber \\[9pt] & = \frac{\Omega_T}{\Omega_{BZ}} \cdot \frac{1} {\varepsilon_{31}\varepsilon_{41}} \left[ \varepsilon_{21}^2 + 3\varepsilon_{21}\varepsilon_{02} + 3\varepsilon_{02}^2 - \frac{\varepsilon_{02}^3\cdot(\varepsilon_{31} + \varepsilon_{42})} {\varepsilon_{32}\varepsilon_{42}} \right] \end{align}
• For $\varepsilon_3 < E \le \varepsilon_4$
$$$n_3(E) = n_2(\varepsilon_3) + \int_{\varepsilon_3}^E \, \rho_3(E)\, \mathrm{d}E \nonumber \\[9pt] = \frac{\Omega_T}{\Omega_{BZ}} \cdot \left[ 1 - \frac{\varepsilon_{40}^3} {\varepsilon_{41}\varepsilon_{42}\varepsilon_{43}} \right]$$$
• For $\varepsilon_4 < E$
$$$n_4(E) = n_3(\varepsilon_4) = \frac{\Omega_T}{\Omega_{BZ}}$$$

### $M_\mathbf{k} \ne 1$

For the case of of general $M_\mathbf{k}$, please refer to the paper by A.H. MacDonald et al., 6 where the deduction was started from $f(\varepsilon_\mathbf{k})$ being a step function and the Dirac delta function case is simply the energy derivative of the former.

### DOS: Smearing vs. Linear-tetrahedron Method

For the calculation of DOS, another often used method is the so-called smearing method, i.e. use a broadening function, for example Gaussian or Lorentzian shape function, to mimic the Dirac function.

As can be seen from Figure 3, with the same density of BZ sampling, the DOS plot from linear-tetrahedron method is much smoother.