Phonon Properties of Diamond-Structure Crystals

for lithium version of Exciting

---

Purpose: In this tutorial, you will learn how to perform a phonon calculation and to obtain properties like frequencies and eigenvectors, the density of modes (or phonon density of states), thermodynamic properties, and phonon-dispersion relations. You will work through the example of diamond.

---

---

0. Prerequisites

This tutorial assumes that you have already downloaded all utility scripts and set the relevant environment variables. Otherwise, please have a look at Tutorial scripts and environment variables. Here, is a list of the scripts which are relevant for this tutorial with a short description.

• THERMAL-cubic.py: Python script for calculating the thermal-expansion coefficient for cubic systems.
• FREE-fromdos.py: Python script for calculating thermodynamic properties.
• PLOT-phonon-dos.gnu: Gnuplot visualization tool for the phonon density of states.
• PLOT-phonon-dispersion.gnu: Gnuplot visualization tool for phonon-dispersion curves.
• PLOT-plot.py: Python visualization tool used by some of the other scripts.

Notes: 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.

In order to start, create a new working directory, e.g., /home/exciting-tutorial/diamond-phonons.

---

1. Set up an input file for a phonon calculation

The setting of the input file can be done in the following steps.

1.1) Define structural parameters

First, the parameter defining your structure should be defined. In this example, we consider carbon in the diamond phase. The structure part of the input file looks then like the following.

  ...
  <structure speciespath="$EXCITINGROOT/species">
 
    <crystal scale="6.7468">
      <basevect>0.5 0.5 0.0</basevect>
      <basevect>0.5 0.0 0.5</basevect>
      <basevect>0.0 0.5 0.5</basevect>
    </crystal>
 
    <species speciesfile="C.xml" rmt="1.25">
      <atom coord="0.00 0.00 0.00" />
      <atom coord="0.25 0.25 0.25" />
    </species>
 
  </structure>
  ...

The chosen lattice parameter of 6.7468 refers to the converged equilibirum parameter for diamond, as obtained with the PBE exchange-correlation functional.

1.2) Parameters for the groundstate calculation

Then, parameters for performing the SCF have to be specified. In this example, we use the GGA functional of Perdew, Burke, and Ernzerhof (PBE), other choices can be found in the description of the attribute xctype of the element groundstate.

  ...
  <groundstate do="fromscratch" 
               ngridk="2 2 2" 
               rgkmax="4.0"
               xctype="GGAPerdew-Burke-Ernzerhof"/>
  ...

Note that this example of input contains values for the parameters ngridk and rgkmax which allow for fast calculations but should not be considered realistic. A detailed description of how to achieve convergence with respect to these parameters can be found in Simple convergence tests.

1.3) Set up the supercell phonon calculation

Simple frozen-phonon calculations can be performed only at high-symmetry points in the Brillouin zone (see, e.g., Phonons at Γ in diamond-structure crystals and Phonons at X in diamond-structure crystals). With exciting, we can also perform phonon calculations at any point in the Brillouin zone using the supercell method and phonon interpolation. The input block in which this method is used is the one specified by the element phonon. In supercell phonon calculations, the size of the supercell is specifyed by the periodicity (the wavevector q) of the phonon displacements which are considered. For instance, as mentioned in Phonons at X in diamond-structure crystals, the supercell corresponding to the X point has a size which is twice the one of the unperturbed crystal.

In general, for calculating phonon properties, you must specify the subelement ngridq inside phonon. Setting ngridq="n1 n2 n3", we are asking to exciting to calculate all the phonons corresponding to all the q-points belonging to a "n1 n2 n3" mesh in reciprocal space. In this case, exciting performs calculation for supercell defined by multiplying the crystal basis vectors of the unit cell by n1, n2, and n3, respectively. The q-points mesh is defined in the same way as the k-points grid (specified by ngridk), which is used in calculating the electronic groundstate. However, the choices of the k-points and q-points meshes are completely uncorrelated.

In the following example, we choose ngridq="2 2 2", which corresponds to a supercell with all crystal basis vectors doubled with respect to the unit cell of the unperturbed crystal.

  ...
  <phonons do="fromscratch" ngridq="2 2 2">
     ...
  </phonons>
  ...

Further attributes of the element phonons are the following.

• deltaph, that controls the step size for atom displacements (default is 0.03 Bohr).
• reduceq, which uses crystal symmetry to reduce the q-points (default is true).

Inside the element phonons, you can insert optional elements, requesting various phonon properties. A basic example is the element qpointset.

  ...
    <qpointset>
        <qpoint> 0.0 0.0 0.0 </qpoint>
        <qpoint> 0.5 0.5 0.0 </qpoint>
        <qpoint> 0.5 0.5 0.5 </qpoint>
    </qpointset>
  ...

This element is needed in order to let exciting calculate phonon eigenvalues (frequencies) and eigenvectors for the q-points specified in each subelement qpoint. In the example above, calculations will be performed for Γ, X, and L phonons (in the standard notation for high-symmetry points in the Brillouin zone of fcc crystals). Note that q-points are given in reciprocal lattice coordinates. Additional properties are discussed in Section 2.

1.4) The basic input file for a phonon calculation

We are now ready to run the actual phonon calculation. Jobs of this type are in general a bit lengthy, the present calculation should take about 30 minutes on a recent CPU with the loose parameters specified. Note that they can also be run in parallel.

<input>
 
  <title>Diamond: Phonons</title>
 
  <structure speciespath="$EXCITINGROOT/species">
 
    <crystal scale="6.7468">
      <basevect>0.5 0.5 0.0</basevect>
      <basevect>0.5 0.0 0.5</basevect>
      <basevect>0.0 0.5 0.5</basevect>
    </crystal>
 
    <species speciesfile="C.xml" rmt="1.25">
      <atom coord="0.00 0.00 0.00" />
      <atom coord="0.25 0.25 0.25" />
    </species>
 
  </structure>
 
  <groundstate do="fromscratch" 
               ngridk="2 2 2" 
               rgkmax="4.0"
               xctype="GGAPerdew-Burke-Ernzerhof"/>
 
  <phonons do="fromscratch" ngridq="2 2 2">
    <qpointset>
        <qpoint> 0.0 0.0 0.0 </qpoint>
        <qpoint> 0.5 0.5 0.0 </qpoint>
        <qpoint> 0.5 0.5 0.5 </qpoint>
    </qpointset>
  </phonons>
 
</input>

Please, remind that after creating the input file input.xml, the path for the species file must be set using the command

$ SETUP-excitingroot.sh

Now, you can run exciting by using the command

$ excitingser

1.5) Output files of a phonon calculation

In addition to the basic output files already mentioned in Electronic-structure calculations, a typical phonon calculation produces the following files.

filename description QPOINTS_PHONON.OUT list of q-points, their coordinates and weights DYN_fileext.OUT complex elements of the dynamical matrix, for every atom (ia) of every species (is) and every polarization direction in crystal coordinates (ip) PHONON.OUT frequencies and eigenvectors of all phonon modes for the q-points specified in the input element qpointset FORCEMAX_fileext.OUT maximum forces, as the phonon input sets automatically tforce to true

Here, fileext means an extension of the filename naming the corresponding q-point of the grid, species, atom, and polarization.

Exercise

• Have a look to the file PHONON.OUT. What is the frequency of the optical phonon modes at Γ?
• What are the eigenvectors of the optical phonon modes at Γ?

---

2. Properties from the phonon calculation

To obtain other phonon properties, additional subelements can be inserted inside the element phonons, as explained in the following.

2.1) Phonon density of states (DOS)

The phonon density of states is computed on top of the phonon calculation by adding the element phonondos within the element phonons.

    ...
    <phonondos ngrdos="100" nwdos="500" ntemp="200" />
    ...

Here, phonon frequencies are computed on a dense grid (ngrdos$\times$ngrdos$\times$ngrdos) and collected to obtain the DOS, which is written to the file PHDOS.OUT. The number of steps between the lowest and highest frequency (N$_{\omega}$) is given by nwdos.

The complete input file in this case is shown below.

<input>
 
  <title>Diamond: PhononDOS</title>
 
  <structure speciespath="$EXCITINGROOT/species">
 
    <crystal scale="6.7468">
      <basevect>0.5 0.5 0.0</basevect>
      <basevect>0.5 0.0 0.5</basevect>
      <basevect>0.0 0.5 0.5</basevect>
    </crystal>
 
    <species speciesfile="C.xml" rmt="1.25">
      <atom coord="0.00 0.00 0.00" />
      <atom coord="0.25 0.25 0.25" />
    </species>
 
  </structure>
 
  <groundstate do="skip" 
               ngridk="2 2 2" 
               rgkmax="4.0"
               xctype="GGAPerdew-Burke-Ernzerhof"/>
 
  <phonons do="skip" ngridq="2 2 2">
    <phonondos ngrdos="100" nwdos="500" ntemp="200" />
  </phonons>
 
</input>

Note the precence of the two attributes do="skip" inside the elements phonons and groundstate, respectively. They let exciting skip a new (and lengthy) phonon calculation and read the necessary data from existing files. This requires you, of course, to have already performed the basic phonon calculation with the same parameters (as in Section 1.4) and to work in the same directory. If this is the case, in order to perform the calculation, you type

$ SETUP-excitingroot.sh
$ excitingser

Data in PHDOS.OUT can be plotted directly. You can use the script PLOT-phonon-dos.gnu, which invokes gnuplot, by typing on the command line

$ PLOT-phonon-dos.gnu

The script produces the PostScript file PLOT.ps displayed here.

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

[show converged plot]

2.2) Thermodynamic properties

The zero-point energy, E_zp, and further thermodynamic properties can be computed from the phonon DOS g$(\omega)$ (see Input reference for further details).

1. the vibrational internal energy E_vib,
2. the vibrational free energy F_vib,
3. the vibrational entropy S_vib
4. the heat capacity c.

In order to calculate and visualize these thermodynamic properties, you can use the script FREE-fromdos.py as follows

$ FREE-fromdos.py 0 1600 40

-----------------------------------
 What do you want to calculate?
-----------------------------------
  0 =>  Zero-point energy
-----------------------------------
  1 =>  Vibrational internal energy
  2 =>  Vibrational free energy
  3 =>  Vibrational entropy
  4 =>  Heat capacity
-----------------------------------

Enter task code >>>> 2

Created PostScript output "PLOT.ps" from file "vibrational-free-energy"

$

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 choice of the thermodynamic quantity that you desire to evaluate (in our example, the vibrational free energy). The result can be visualized by watching the PostScript output file PLOT.ps.

Note that if the calculation of the zero-point energy is required, this energy is not plotted, but it is written in output directly on the screen, e.g.,

...
Zero-point energy is  1.37898903e-02  [Ha]
...

The converged plot for the vibrational free energy can be seen here.

[show converged plot]

Note the difference in the energy scale of the converged plot.

2.3) Phonon-dispersion relations

Phonon-dispersion relations can also be easily calculated. Phonon eigenvalues are computed along a path through the Brillouin zone and written for later plotting. This path can be specified within the element phonondispplot (reciprocal lattice coordinates) in the following way (a path between the points Γ-(K)-X-Γ-L is specified in the present example).

    ...
    <phonondispplot>
      <plot1d>
        <path steps="200">
          <point coord="1.0   0.0   0.0" />
          <point coord="0.5   0.5   0.0" />
          <point coord="0.0   0.0   0.0" />
          <point coord="0.5   0.0   0.0" />
        </path>
      </plot1d>
    </phonondispplot>
    ...

This codelet is to be included within the element phonons, same as for qpointset and phonondos above. Of course, one can reuse data of a previous phonon calculation by adding do="skip" to the elements phonons and groundstate. This run produces the file PHDISP.OUT that contains the phonon eigenvalues along the path specified in the input and can be plotted directly. The script PLOT-phonon-dispersion.gnu can be used to produce a plot of the phonon-dispersion relations.

$ PLOT-phonon-dispersion.gnu 0 1600

Here, you must also specify the minimum and maximum frequency in the plot. The resulting PostScript file PLOT.ps looks like the following.

Converged results as mentioned above can be seen here.

[show converged plot]

2.4) Other useful features

[show more features]

• The element interpolate triggers an interpolation on a given q-point grid. For every point in the grid, which may be offset from the origin, phonon frequencies and (optionally) eigenvectors can be obtained, and written to the file PHONON_INTERPOLATE.OUT.
• The element reformatdynmat (without any attributes) requests an output of the dynamical matrix in the form of a 3$\times$3 block for each combination of atoms and each q-point (file DYNMAT.OUT).
• The file DYNMAT_SUMRULE.OUT contains the dynamical matrix (in 3x3 blocks as above) after applying a correction, such that the matrix for q=0 satisfies the acoustic sum rule. This ensures that the q=0 dynamical matrix has 3 zero eigenvalues. The files DYNMAT_SYM.OUT and DYNMAT_SYM_SUMRULE.OUT contain the symmetrized matrices.

---

3. Thermal expansion

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 show how to calculate with exciting physical quantities which are related to the phenomenon of thermal expansion.

3.1) Mode Grüneisen parameters

Phonon frequencies of a real crystal explicitly depend on the crystal volume. 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

(1)

\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

(2)

\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.

In order to calculate the mode Grüneisen parameter at some high-symmetry point in the first Brillouin zone of the diamond lattice, create two new working directories, e.g., diamond-phonons-minus and diamond-phonons-plus. Copy the input file already used for calculating the phonon dispersion at the (zero-temperature) equilibrium lattice parameter $a_0$ to each of the two directories. In diamond-phonons-minus change the lattice parameter in the input file to $a_0-\Delta$, with $\Delta$=0.05 Bohr. Similarly, change the lattice constant in the input file in diamond-phonons-plus to $a_0+\Delta$.

Now, repeat the calculation of the phonon frequencies of diamond at Γ and X as above in each of the two new directories. Compute the mode Grüneisen parameters using the above definition (evaluating the derivative numerically from a quadratic interpolation, i.e., from the differences between $\omega_j(a_0 + \Delta)$ and $\omega_j(a_0 - \Delta)$). 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}$ 1611.76 1574.45 1537.86 1.06 0.99 0.96 $X_{\rm TA}$ 721.51 722.72 723.07 -0.05 0.22 $X_{\rm TO}$ 1177.33 1133.34 1089.89 1.74 1.51 $X_{\rm L}$ 1180.08 1152.70 1125.89 1.06 0.91

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

3.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

(3)

\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

(4)

\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 with respect to the volume. The derivative is approximated numerically by a quadratic interpolation.

(5)

\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 Section 3.1.
2. Create the new directory thermal-expansion.
3. Copy the files input.xml and thermo.xml from the directory diamond-phonon 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 6$\times$6$\times$6 k-grid, a rgkmax of 7.0 and 4 q-points in each direction) can be found here.

[show converged plot]

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

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.

---