by Pasquale Pavone & Stefan Kontur for exciting beryllium
Purpose: In this tutorial, you will learn how to perform a series of phonon calculations in order to obtain the mode Grüneissen 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 Tutorial scripts and environment variables. Furthermore, we also assume that you are familiar with the concepts and the calculations presented in Phonon properties of diamondstructure crystals. The only scripts which is relevant for this tutorial is
 THERMALcubic.py: Python script for calculating the thermalexpansion 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^{5} 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)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)where $V_0$ is the equilibrium volume. For cubic systems and in terms of the lattice parameter $a$ one can write
(3)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 thermalexpansion coefficient $\alpha(T)$. Within the quasiharmonic approximation, one can write
(4)where $B_0$ is the bulk modulus and
(5)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 highsymmetry 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. (25). 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.,
 diamondphononsminus,
 diamondphononszero,
 diamondphononsplus,
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 diamondstructure crystals, you can reduce the computational effort by copying the results you have obtained in the "equilibrium" directory diamondphononszero. If this is not the case, create an input file for the phonon calculation as shown in Phonon properties of diamondstructure crystals inside each of the three directories. Change the lattice parameter in the input file inside diamondphononsminus to $a_00.05$. Similarly, change the lattice constant in the input file in diamondphononsplus 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 diamondstructure 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)The following values should be obtained (frequencies are given in cm^{1}).

Here, $\gamma_{\rm conv}$ is the converged result, obtained using a 8$\times$8$\times$8 kpoint mesh, rgkmax = "7.0", and a grid of 4$\times$4$\times$4 qpoints.
2.2) Thermal expansion coefficient
The derivative in Eq. (5) is approximated numerically by a quadratic interpolation.
(7)In order to calculate $\alpha(T)$, you can use the following procedure.
 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 diamondstructure crystals.
 Create the new directory thermalexpansion.
 Copy the files input.xml and thermo.xml from the directory diamondphononzero to the directory thermalexpansion adding the suffix .0. Thus, you will have inside thermalexpansion the two new files input.xml.0 and thermo.xml.0.
 Repeate the previous step from the directory diamondphononplus and the suffix .+.
 Repeate the previous step from the directory diamondphononminus and the suffix ..
Now, inside the directory thermalexpansion you have the following files.
input.xml. thermo.xml.
input.xml.0 thermo.xml.0
input.xml.+ thermo.xml.+
The linear thermalexpansion coefficient can be obtained using the script THERMALcubic.py inside the directory thermalexpansion.
$ THERMALcubic.py 0 1600 40

Linear thermalexpansion 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 kgrid, rgkmax = "7.0", and 4 qpoints in each direction) can be found below.
Exercise: Silicon in diamond structure
 The calculation presented here for diamond can be repeated for Silicon. Use the PBE exchangecorrelation functional, an equilibrium lattice parameter of 10.34 Bohr, and low values for both k and qpoint grids to accelerate the runs, as above. The experimental bulk modulus is 98 GPa (LandoltBö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.