Phonon and Thermal Properties in Diamond-Structure Systems

Purpose: In this tutorial you will learn how to calculate phonon-dispersion relations and some thermal properties in crystals with the diamond structure. As an example we consider Si in the diamond structure.

1. Define relevant shell variable and download scripts

Very important: Before starting, the following shell variable must be set by the user:

EXCITINGROOT = Directory where exciting has been downloaded, e.g.: /home/user/exciting .

The setting of this variable can be done in a bash shell by typing (from now on the symbol $ will indicate the shell prompt):

$ export EXCITINGROOT=/local_path_to_exciting_root

The directory name /local_path_to_exciting_root is a dummy name which must be explicitly changed by the user to the appropriate one.

Please note: In the input-files shown on this page, the placeholder EXCITINGROOT always needs to be replaced by the value of the $EXCITINGROOT variable.

As a second step, it is necessary to download the scripts which are used in this tutorial. Here is a list of these scripts with a short description, click on the script name to download it:

  • GNU-dos: Gnuplot visualization tool for the phonon density of states.
  • GNU-phdisp: Gnuplot visualization tool for the phonon-dispersion relations at a given lattice parameter.
  • GNU-3phdisp: Gnuplot visualization tool for the phonon-dispersion relations at at three different lattice parameters.
  • GNU-thermo: Gnuplot visualization tool for the fit of the vibrational free and internal energy as a function of temperature.

2. Set up the input file

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

i) Define structural parameters

First, the parameter defining your structure should be defined. In this example, we consider the case of Si in the diamond structure, then the structure part of the input file should look like the following:

  ...
  <structure speciespath="EXCITINGROOT/species">
 
    <crystal scale="10.21">
      <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="Si.xml">
      <atom coord="0.00 0.00 0.00" />
      <atom coord="0.25 0.25 0.25" />
    </species>
 
  </structure>
  ...
ii) Ground-state calculation

Then, the parameters for performing the SCF should be specified as in the example:

  ...
  <groundstate do="fromscratch" 
               ngridk="2 2 2" 
               findlinentype="simple"
               lmaxapw="8"
               epsband="1.0d-5" 
               lmaxvr="7"/>
  ...
iii) Set super-cell phonon calculations

In order to perform an exciting calculation which can takes into account the periodicity of a phonon with wave-vector q, one has to choose an appropriate supercell. For instance, as mentioned in the tutorial 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 for all the q-points belonging to a given mesh in the Brillouin zone, say the one corresponding to ngridq="n1 n2 n3", exciting will perform calculation for supercell defined by multiplying the crystal basis vectors of the unit cell by n1, n2, and n3, respectively.

In the following example, we choose the value ngridq="2 2 2", which correspond 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>
  ...
iv) Set calculation of the phonon density of states

Similarly to the electronic density of states (see the tutorial Electronic-structure calculations) the calculation of the phonon DOS is specified by:

    ...
    <phonondos/>
    ...
v) Set calculation of the phonon dispersion relations
    ...
    <phonondispplot>
      <plot1d>
        <path steps="200">
          <point coord="0.0   0.0  -1.0" label="GAMMA" />
          <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.5   0.5" label="L"     />
        </path>
      </plot1d>
    </phonondispplot>
    ...
vi) The complete input file
<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?>
 
<input xsi:noNamespaceSchemaLocation="EXCITINGROOT/excitinginput.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" 
  xsltpath="EXCITINGROOT/xml/">
 
  <title>Si: Phonon</title>
 
  <structure speciespath="EXCITINGROOT/species">
 
    <crystal scale="10.21">
      <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="Si.xml">
      <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" 
               findlinentype="simple"
               lmaxapw="8"
               epsband="1.0d-5" 
               lmaxvr="7"/>
 
  <phonons do="fromscratch" ngridq="2 2 2">    
 
    <phonondos/>
 
    <phonondispplot>
      <plot1d>
        <path steps="200">
          <point coord="0.0   0.0  -1.0" label="GAMMA" />
          <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.5   0.5" label="L"     />
        </path>
      </plot1d>
    </phonondispplot>
 
  </phonons>
 
</input>
si-dos-q222-a10.21.png
si-phdisp-q222-a10.21.png
si-dos-q444-a10.21.png
si-phdisp-q444-a10.21.png 
  ...
  <groundstate do="skip" 
    ...     
  <phonons do="skip" ngridq="2 2 2">    
    ...

<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?>
 
<input xsi:noNamespaceSchemaLocation="EXCITINGROOT/excitinginput.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" 
  xsltpath="EXCITINGROOT/xml/">
 
  <title>Si: Phonon</title>
 
  <structure speciespath="EXCITINGROOT/species">
 
    <crystal scale="10.21">
      <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="Si.xml">
      <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" 
               findlinentype="simple"
               lmaxapw="8"
               epsband="1.0d-5" 
               lmaxvr="7"/>
 
  <phonons do="skip" ngridq="2 2 2">    
 
    <phonondos/>
 
    <phonondispplot>
      <plot1d>
        <path steps="200">
          <point coord="0.0   0.0  -1.0" label="GAMMA" />
          <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.5   0.5" label="L"     />
        </path>
      </plot1d>
    </phonondispplot>
 
  </phonons>
 
</input>
si-3disp-q444.png

3. Thermodynamic properties

The thermodynamic properties discussed here are all obtained from the phonon density of states (DOS). The following quantities are calculated:

  • Zero point energy
(1)
\begin{align} E_{\rm ZP} = \frac{1}{2} N_{\rm at} \int_{\omega_{\rm min}}^{\omega_{\rm max}} \omega g(\omega) \end{align}
  • Vibrational energy
(2)
\begin{align} E_{\rm vib} = \frac{1}{2} N_{\rm at} \int_{\omega_{\rm min}}^{\omega_{\rm max}} \omega g(\omega) \coth\frac{\omega}{2k_B T} \end{align}
  • Vibrational free energy
(3)
\begin{align} F_{\rm vib} = N_{\rm at} k_B T \int_{\omega_{\rm min}}^{\omega_{\rm max}} g(\omega) \log\left(2 \sinh \frac{\omega}{2k_B T}\right) \end{align}
  • Vibrational entropy
(4)
\begin{align} S_{\rm vib} = \frac{E_{\rm vib}-F_{\rm vib}}{T} \end{align}
  • Heat capacity
(5)
\begin{align} c = N_{\rm at} k_B \int_{\omega_{\rm min}}^{\omega_{\rm max}} g(\omega) \left(\frac{\omega}{k_B T}\right)^2 \frac{\exp(\frac{\omega}{k_B T})}{\left[\exp(\frac{\omega}{k_B T})-1\right]^2} \end{align}

For these expressions above the following abbreviations are used

  • Phonon DOS
(6)
\begin{align} g(\omega) \end{align}
  • Temperature grid
(7)
\begin{align} T_{\rm max} = \omega_{\rm max} / k_B \end{align}
(8)
\begin{align} T_j = T_{\rm max} \frac{j}{N_T} \end{align}
  • $N_T$ is the number of Temperatures for thermodynamical properties.
  • The range of phonon frequencies is defined as follows
(9)
\begin{align} \omega_{\rm max} = {\rm max}_j\{\omega_j\} \end{align}
(10)
\begin{align} \omega_{\rm min} = {\rm min}_j\{\omega_j\} \end{align}
(11)
\begin{align} \Delta\omega = \omega_{\rm max} - \omega_{\rm min} \end{align}
(12)
\begin{align} \delta\omega = \Delta\omega / N_{\omega} \end{align}
  • Here $N_{\omega}$ is the number of frequency points for sampling the phonon DOS.
  • The range of phonon frequencies is slightly enlarged
(13)
\begin{align} \omega_{\rm max} \rightarrow \omega_{\rm max} +\frac{\delta \omega}{N_T} \end{align}
(14)
\begin{align} \omega_{\rm min} \rightarrow \omega_{\rm min} -\frac{\delta \omega}{N_T} \end{align}

Results

si-thermo-q444-a10.21.pngsi-cv-q444-a10.21.png

The vibrational entropy is obtained by issuing the following command

xsltproc $EXCITINGROOT/xml/visualizationtemplates/xmlthermo2agr.xsl thermo.xml

which will generate the following grace-files

heat_capacity.agr
vibrational_energy.agr
vibrational_entropy.agr
vibrational_free_energy.agr

The vibrational entropy for a=10.26

si-vibentrop-q444-k666-a10.26.png
exercisetutorial
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License