Thermal Expansion of Diamond-Structure Crystals

by Pasquale Pavone & Stefan Kontur for exciting boron

Purpose: In this tutorial, you will learn how to perform a series of phonon calculations in order to obtain the mode Grüneisen parameters and the thermal-expansion coefficient of diamond.

---

0. Prerequisites

This tutorial assumes that you have already set the relevant environment variables. Otherwise, please have a look at How to set environment variables for tutorials scripts. Furthermore, we also assume that you are familiar with the concepts and the calculations presented in Phonon properties of diamond-structure crystals. The only scripts which is relevant for this tutorial is

• THERMAL-cubic.py: Python script for calculating the thermal-expansion coefficient for cubic systems.

Note: The symbol $ at beginnings of lines in code segments below indicates the shell prompt.

In the following we use the conversion factor 1 Ha $\approx$ 2.194 746 x 10⁵ cm^-1.

All results obtained in this tutorial for diamond can be compared with theoretical and experimental data present in the literature.

---

1. Theoretical background

Harmonic phonon frequencies do not depend on the crystal volume. As a consequence of this fact, there is no effect of temperature on the equilibrium volume of a perfect harmonic crystal. However, real crystals do change their volume as their temperature is varied. In the next, we introduce physical quantities which are related to the phenomenon of thermal expansion.

1.1) Mode Grüneisen parameters

Phonon frequencies of a real crystal explicitly depend on the crystal volume

(1)

\begin{align} \omega_j({\bf q}) = \omega_j({\bf q},V)\,. \end{align}

This dependence can be characterized by the mode Grüneisen parameters, $\gamma_j({\bf q})$. The Grüneisen parameter of the mode $j$ at a wavevector ${\bf q}$ is defined as

(2)

\begin{align} \gamma_j({\bf q}) = -\frac{V_0}{\omega_j({\bf q},V_0)}\left(\frac{\partial\,\,\omega_j({\bf q},V)}{\partial\, V}\right)_{V=V_0}\,, \end{align}

where $V_0$ is the equilibrium volume. For cubic systems and in terms of the lattice parameter $a$ one can write

(3)

\begin{align} \gamma_j({\bf q}) = -\frac{a_0}{3\,\,\omega_j({\bf q},a_0)}\left(\frac{\partial\,\,\omega_j({\bf q},a)}{\partial\, a}\right)_{a=a_0}\,, \end{align}

where again $a_0$ is the equilibrium parameter.

1.2) Thermal expansion coefficient

The thermal expansion of a crystal can be quantified by the linear thermal-expansion coefficient $\alpha(T)$. Within the quasi-harmonic approximation, one can write

(4)

\begin{align} \alpha(T) = \frac{1}{3V}\frac{\partial V}{\partial T} \approx -\frac{1}{3V_0B_0}\frac{\partial\,\varphi(T)}{\partial T}\,, \end{align}

where $B_0$ is the bulk modulus and

(5)

\begin{align} \varphi(T) = \left(\frac{\partial F_{\rm vib}(V,T)}{\partial V}\right)_{V=V_0} \end{align}

is the derivative of the vibrational free energy $F_{\rm vib}(V,T)$ with respect to the volume.

---

2. Calculations

The calculation of the thermal expansion coefficient and of the mode Grüneisen parameter at some high-symmetry point in the first Brillouin zone of the diamond lattice requires the evaluation of derivatives with respect to the volume or, equivalently, to the lattice constant, as shown in Eqs. (2-5). In order to calculate these derivatives in a finite difference approach, phonon calculations have to be performed at three different volumes (lattice constants).

Therefore, the first step is to create three new working directories, e.g.,

• diamond-phonons-minus,
• diamond-phonons-zero,
• diamond-phonons-plus,

which will correspond to the calculations performed at lattice parameters $a_0-\Delta$, $a_0$, and $a_0+\Delta$, respectively. Here, $a_0$ is the equilibrium lattice parameter and $\Delta$ a positive "small" number which in this tutorial is set to 0.05 Bohr.

If you have already performed the calculation presented in Phonon properties of diamond-structure crystals, you can reduce the computational effort by copying the results you have obtained in the "equilibrium" directory diamond-phonons-zero. If this is not the case, create an input file for the phonon calculation as shown in Phonon properties of diamond-structure crystals inside each of the three directories. Change the lattice parameter in the input file inside diamond-phonons-minus to $a_0-0.05$. Similarly, change the lattice constant in the input file in diamond-phonons-plus to $a_0+0.05$.

2.1) Mode Grüneisen parameters

Now, repeat the calculation of the phonon frequencies of diamond at Γ and X as mentioned in Section 1.4 of Phonon properties of diamond-structure crystals inside each new directory. When all executions are completed, compute the mode Grüneisen parameters using Eq. (3) and evaluating the numerical derivative from a quadratic interpolation,

(6)

\begin{align} \left.\frac{\partial\,\,\omega_j({\bf q},a)}{\partial\, a}\right|_{a=a_0}\approx \frac{\omega_j({\bf q},a+\Delta) - \omega_j({\bf q},a-\Delta)}{2\Delta}\,. \end{align}

The following values should be obtained (frequencies are given in cm^-1).

q-point $\omega(a_0 - \Delta)$ $\quad\omega(a_0)$ $\omega(a_0 + \Delta)$ $\quad\gamma$ $\gamma_{\rm conv}$ $\gamma_{\rm expt}$ $\Gamma_{\rm LTO}$ 1618.43 1580.69 1543.78 1.06 1.00 0.96 $X_{\rm TA}$ 511.86 527.51 541.26 -1.25 0.21 $X_{\rm TO}$ 1126.14 1097.96 1070.56 1.14 1.52 $X_{\rm L}$ 1174.86 1128.90 1083.77 1.81 0.91

Here, $\gamma_{\rm conv}$ is the converged result, obtained using a 8$\times$8$\times$8 k-point mesh, rgkmax = "7.0", and a grid of 4$\times$4$\times$4 q-points.

2.2) Thermal expansion coefficient

The derivative in Eq. (5) is approximated numerically by a quadratic interpolation.

(7)

\begin{align} \varphi(T) \approx \frac{2}{3} \frac{F_{\rm vib}(a_0+\Delta,T) - F_{\rm vib}(a_0-\Delta,T)}{a_0^2\,\Delta} \end{align}

In order to calculate $\alpha(T)$, you can use the following procedure.

1. Obtain the phonon DOS and thermodynamical properties for $a_0 - \Delta$, $a_0$, and $a_0 + \Delta$ using the same procedure of Sections 2.1 and 2.2 of Phonon properties of diamond-structure crystals.
2. Create the new directory thermal-expansion.
3. Copy the files input.xml and thermo.xml from the directory diamond-phonon-zero to the directory thermal-expansion adding the suffix .0. Thus, you will have inside thermal-expansion the two new files input.xml.0 and thermo.xml.0.
4. Repeate the previous step from the directory diamond-phonon-plus and the suffix .+.
5. Repeate the previous step from the directory diamond-phonon-minus and the suffix .-.

Now, inside the directory thermal-expansion you have the following files.

input.xml.-     thermo.xml.-
input.xml.0     thermo.xml.0
input.xml.+     thermo.xml.+

The linear thermal-expansion coefficient can be obtained using the script THERMAL-cubic.py inside the directory thermal-expansion.

$ THERMAL-cubic.py 0 1600 40

------------------------------------------------------------------------
Linear thermal-expansion coefficient for cubic systems
------------------------------------------------------------------------

Enter value for the bulk modulus [GPa] >>>> 452

$

where the 3 online entries are the minimum temperature (usually 0), the maximum temperature, and the number of temperature steps between the minimum and the maximum. The script also requires the value of the bulk modulus, in our example the experimental value for diamond has been used (452 GPa). Alternatively, the bulk modulus can be obtained as described in Volume optimization for cubic systems. The result can be visualized by watching the PostScript output file PLOT.ps.

A converged result (using an 8$\times$8$\times$8 k-grid, rgkmax = "7.0", and 4 q-points in each direction) can be found below.

[show converged plot]

Exercise: Silicon in diamond structure

• The calculation presented here for diamond can be repeated for Silicon. Use the PBE exchange-correlation functional, an equilibrium lattice parameter of 10.34 Bohr, and low values for both k- and q-point grids to accelerate the runs, as above. The experimental bulk modulus is 98 GPa (Landolt-Börnstein tables).
• What are the differences to diamond?
• Use the theoretical bulk modulus. Calculate theoretical values of equilibrium parameters as in Volume Optimization for Cubic Systems.
• Compare your results with the "converged" one displayed here.