Lattice Dynamics of Diamond and Zincblende-Structure Crystals

by Sebastian Tillack & Pasquale Pavone for exciting neon

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 examples of diamond and cubic Boron nitride.

---

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. Here, is a list of the scripts which are relevant for this tutorial with a short description.

• FREE-fromdos.py: Python script for calculating thermodynamic properties.
• PLOT-dos.py: Python visualization tool for the phonon density of states.
• PLOT-band-structure.py: Python visualization tool for phonon-dispersion curves.
• PLOT-files.py: Python visualization tool used by some of the other scripts.

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.

---

1. Theoretical Background

---

☛ Theory of linear response for phonons

Within density-functional perturbation theory (DFPT) we calculate the linear response of the Kohn-Sham (KS) system upon an external perturbation $\delta$. We express the quantities of the perturbed system as an expansion in the perturbation operator and write down the KS equations for the perturbed system

(1)

\begin{align} \left(\hat{\bf h}\phantom{}^0 + \lambda \,\delta \hat{\bf h} + \dots\right) \left(\psi_{n{\bf k}}^0({\bf r}) + \lambda\, \delta \psi_{n{\bf k}}({\bf r}) + \dots\right) = \left(\epsilon_{n{\bf k}}^0 + \lambda \,\delta \epsilon_{n{\bf k}} + \dots\right) \left(\psi_{n{\bf k}}^0({\bf r}) + \lambda \,\delta \psi_{n{\bf k}}({\bf r}) + \dots\right) \;. \end{align}

Collecting terms linear in $\lambda$ we find the so called Sternheimer equation

(2)

\begin{align} \left[\hat{\bf h}\phantom{}^0 - \epsilon_{n{\bf k}}^0\right] \,\delta \psi_{n{\bf k}}({\bf r}) = -\left[\delta \hat{\bf h} - \delta \epsilon_{n{\bf k}}\right] \,\psi_{n{\bf k}}^0({\bf r}) \;, \end{align}

from which the linear response of the wavefunctions $\delta \psi_{n{\bf k}}$ can be obtained. Equivalently to the KS equations, the Sternheimer equation has to be solved self consistently, since the Hamiltonian response implicitly depends on the wave function response via the density and potential response:

(3)

\begin{align} \delta \psi_{n{\bf k}} \rightarrow \delta n({\bf r}) \rightarrow \delta V({\bf r}) \rightarrow \delta \hat{\bf h} \rightarrow \delta \psi_{n{\bf k}} \;. \end{align}

In the case of DFPT for phonons, the external perturbation is a collective phonon-like displacement of the nuclei from their equilibrium positions $\pmb{\tau}_{\kappa{\bf R}}$

(4)

\begin{align} \delta = \delta^{\bf q}_{\kappa \alpha} \quad:\quad \pmb{\tau}_{\kappa {\bf R}} \longrightarrow \pmb{\tau}_{\kappa {\bf R}} + {\rm e}^{{\rm i}{\bf q}\cdot{\bf R}} \delta \pmb{\tau}_{\kappa} \;, \end{align}

where $\kappa$ and ${\bf R}$ label atoms and unit-cells, respectively, ${\bf q}$ is the phonon wave vector and $\delta \pmb{\tau}_{\kappa}$ is an infinitesimal displacement vector.

In order to exploit symmetries, we do not displace every atom individually in the three Cartesian directions (canonical basis) but we use special symmetry adapted displacement patterns (irreducible representations, irreps):

(5)

\begin{align} \delta^{\bf q}_{Id} = \sum_{\kappa,\alpha} p_{\kappa \alpha, Id}({\bf q})\, \delta^{\bf q}_{\kappa \alpha} \;, \end{align}

where $I$ labels an irrep (for a given ${\bf q}$) and $d$ labels a so-called member of the irrep.

Once we have converged the Sternheimer equation for the density and potential response, we can use the results to compute in a single step the linear response of the atomic forces ${\bf F}_{\kappa}$ which directly correspond to the entries of the dynamical matrix

(6)

\begin{align} \delta^{\bf q}_{\kappa \alpha} F_{\kappa' \beta} = \sqrt{M_\kappa\, M_{\kappa'}}\, D_{\kappa \alpha, \kappa' \beta}({\bf q}) \;. \end{align}

---

2. Phonon properties for diamond

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

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

$ cd /home/tutorials
$ mkdir diamond-phonons
$ cd diamond-phonons

---

2.1 Set-up the input file and run calculation

Create the input.xml file with the following content.

<input>
 
   <title>Diamond: Phonons</title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="6.749908">
         <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="4 4 4"
      swidth="0.0001"
      rgkmax="4.0"
      gmaxvr="14"
      tforce="true">
      <libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
   </groundstate>
 
   <phonons
      do="fromscratch"
      method="dfpt"
      ngridq="2 2 2"
      polar="false">
      <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>

The input block that describes a phonon calculation is specified by the element phonons. In exciting, there are two different methods available to compute phonons: the supercell method and DFPT. The latter is the default and what will cover in the first part of this tutorial. The dynamical matrices are calculated for a regular grid of phonon wavevectors q specified by the attribute ngridq inside phonons. 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 ground state. Afterwards, the phonon frequencies and eigenvectors at any wavevector q can be obtained by Fourier interpolation. Note that in the case of a DFPT calculation, the k- and q-grid need to be commensurate, i.e., ngridk must be an integer multiple of ngridq. In the case of a supercell calculation, both grids are completely uncorrelated.

A novelty in the groundstate element is that we replace the attribute xctype by the new element libxc. By that, we specify to take the exchange-correlation potential from the libxc library instead of the exciting internal implementation. The reason for this change is that a DFPT calculation requires the xc-kernel (the second derivatives of the xc-energy) which is not natively implemented in exciting and is always obtained from the libxc library. For full consistency, we recommend to also use the xc-potential (first derivative of the xc-energy) from libxc. However, in principle you can also use the exciting internal implementations of the xc-potential using the attribute xctype. In this case, only the xc-kernel of the corresponding potential is obtained from libxc.

Further attributes of the element phonons are the following.

• method, that controls the method employed to compute phonons (default is "dfpt").
• reduceq, which uses crystal symmetry to reduce the q-points (default is "true").
• polar, that specifies whether a polar or nonpolar material is to be calculated (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.3.

After you have created the input file input.xml, you can run exciting by using the commands

$ SETUP-excitingroot.sh
$ time exciting_smp &

The calculation is completed if the file PHONON.OUT has been created. This can be checked by the command

$ ls PHONON.OUT

---

2.2 Output files of a phonon calculation

The main result of a phonon calculation are the dynamical matrices. They are given in the files DYN_fileext.OUT, where fileext means an extension of the filename having the structure Q####_####_####_S##_A###_P#.OUT and specifying the q vector on the grid (Q…) species and index of the displaced atom (S… and A…) and polarization direction (P…). Additionally, a number of subdirectories is created. In a DFPT calculation, their names are Q####_####_####_I###, where the first part specifies the q-point in the grid and the second part labels the irrep at that specific q-point. Irreps are special symmetry adapted displacement patterns. Inside of each subdirectory, you will find the following files:

filename description PHONON_INFO_fileext.OUT General information about the run. Similar to INFO.OUT for the groundstate. DYN_fileext_#.OUT Force response (or dynamical matrix row) for the #-th member of the irrep. PHONON_DRHO_fileext_#.OUT Electron density response for the #-th member of the irrep. PHONON_DVEFF_fileext_#.OUT Effective potential response for the #-th member of the irrep.

where fileext is the same as the subdirectory name.

In the following, we take a closer look to the file PHONON_INFO_Q0000_0000_0000_I001.OUT.

2.2.1 PHONON_INFO_fileext.OUT

The first part of the file contains information about the use of symmetries, parallelization settings and the whole phonon run.

********************************************************************************
* symmetries and parallelization settings                                      *
********************************************************************************

Irreducible representations are used.
k-points are reduced by symmetries.
Response quantities are symmetrized.

q-points and symmetries:

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    8        2 (3*,3)
   2 ( 0.500000, 0.000000, 0.000000)     12   13        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16   13        3 (2,2,2)

current task / total (remaining) number of tasks:         1/       9 (9)
load of current task / total (remaining) load:           24/     204 (204)
procs for current task / total number of procs:           1/       1
average / maximum load per proc:                        204/     204
number of idle procs:                                     0

The first lines tell that the irreps have been used as displacement patterns and that the computational cost has been reduced by exploiting symmetries.

The following table displays all (irreducible) q-points that are computed in that phonon calculation. N_sym is the number of crystal symmetries in the small group of q. N_k indicates the number of k-points in the irreducible Brillouin zone obtained with the symmetries in the small group. N_irrep denotes the number of irreps for the given q-point and the values in parenthesis give their respective size (or equivalently the number of members per irrep). The q-point and irrep this file corresponds to is marked with an asterisk.

From the last lines we can extract that the whole calculation consists of 9 independent tasks from which all 9 are remaining and that this particular file corresponds to task 1. The meaning of the other lines will be discussed below in the Section 5.2.

In the next part, the displacement patterns of the irreps are given.

********************************************************************************
* irreducible representations / atomic displacement patterns                   *
********************************************************************************

irrep member 1
  C     1    +0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
  C     2    -0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i

irrep member 2
  C     1    +0.0000000000 +0.0000000000i
             +0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
  C     2    +0.0000000000 +0.0000000000i
             -0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i

irrep member 3
  C     1    +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.7071067812 +0.0000000000i
  C     2    +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
             -0.7071067812 +0.0000000000i

In this example, the irrep has three members and for each member, the three components of the (potentially complex) displacement vector for each atom are given.

The remainder of the file contains information about the self-consistency cycle of the Sternheimer equation and shows the change in the potential response for each iteration:

********************************************************************************
* scf loop for Sternheimer equation started                                    *
********************************************************************************

Sternheimer SCF loop iteration   1
maximum change in potential response:                0.960508E-01
Fermi energy response:                                0.00000    
Time spent for iteration (seconds):                          1.20

Sternheimer SCF loop iteration   2
maximum change in potential response:                0.646932E-01
Fermi energy response:                                0.00000    
Time spent for iteration (seconds):                          3.40

...

2.2.2 PHONON.OUT

If you have specified a qpointset in the input file, the frequencies of the phonon modes at the required q-points are listed inside the file PHONON.OUT. However, in this file frequencies are given in non frequently used atomic units. In order to get frequency values which can be compared to experimental data, you may use the script CONVERT-au2invcm.py to obtain frequencies expressed in units of cm^-1.

$ CONVERT-au2invcm.py

############################################

--------------------------------------------
Total number of atoms    :     2
Total number of q-points :     3
--------------------------------------------

q-point   1  is  0.0 0.0 0.0
   mode   1      frequency:    -0.0000  cm-1
   mode   2      frequency:    -0.0000  cm-1
   mode   3      frequency:    -0.0000  cm-1
   mode   4      frequency:  1244.0005  cm-1
   mode   5      frequency:  1244.0005  cm-1
   mode   6      frequency:  1244.0005  cm-1

q-point   2  is  0.5 0.5 0.0
   mode   1      frequency:   782.9265  cm-1
   mode   2      frequency:   782.9265  cm-1
   mode   3      frequency:  1006.6311  cm-1
   mode   4      frequency:  1006.6311  cm-1
   mode   5      frequency:  1195.6297  cm-1
   mode   6      frequency:  1195.6297  cm-1

q-point   3  is  0.5 0.5 0.5
   mode   1      frequency:   539.9045  cm-1
   mode   2      frequency:   539.9045  cm-1
   mode   3      frequency:   998.7204  cm-1
   mode   4      frequency:  1149.7671  cm-1
   mode   5      frequency:  1149.7671  cm-1
   mode   6      frequency:  1275.9419  cm-1

############################################

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.3 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.3.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 ngridqint="70 70 70" nwdos="1000" nsmdos="2"/>
...

Here, phonon frequencies are computed on a dense grid (ngridqint) 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, nsmdos is the number of times the data are smoothened.

The modified input file in this case is shown below.

<input>
 
   <title>Diamond: Phonons</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      ... 
   </structure>
 
   <groundstate 
      do="skip" 
      .../>
 
   <phonons 
      do="skip"
      ...>
      <phonondos 
         ngridqint="70 70 70" 
         nwdos="1000"  
         nsmdos="2"/>
   </phonons>
 
</input>

Note the presence 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 2.1) and to work in the same directory. If this is the case, in order to perform the calculation, you type

$ SETUP-excitingroot.sh
$ time exciting_smp &

Data in PHDOS.OUT can be plotted directly. You can use the script PLOT-dos.py (see The python script "PLOT-dos.py"), by typing on the command line

$ PLOT-dos.py -p  -e 0 1400

The script produces a PNG (PLOT.png) output, which should be similar to the one displayed here.

A better 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

Hint: Run converged calculation in a different directory, e.g., /home/tutorial/diamond-phonons/conv.

$ mkdir /home/tutorial/diamond-phonons/conv
$ cd /home/tutorial/diamond-phonons/conv

Copy here the previous input.xml and modify the computational parameters. The result should look like the following.

---

2.3.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 free energy F_vib= U_vib -TS_vib,
2. the vibrational internal energy U_vib,
3. the entropic contribution to the vibrational free energy TS_vib,
4. the vibrational entropy S_vib
5. the heat capacity C_v.

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

$ FREE-fromdos.py 0 1400 160 eV

 Zero-point energy is 3.5124e-01 [eV]

 Created PNG output "PLOT-free-energies.png"
 Created PNG output "PLOT-entropy.png"
 Created PNG output "PLOT-heat-capacity.png"

where the 4 online entries are the minimum temperature (usually 0), the maximum temperature, the number of temperature steps between the minimum and the maximum, and the energy units (either eV or Ha) . The value of the calculated zero-point energy is written on the screen output. The script also generates three PNG output files: PLOT-free-energies.png, PLOT-entropy.png, and PLOT-heat-capacity.png. These results can be visualized by watching directly the PNG files with standard tools. Examples are given in the following:

Better converged plots can be seen here.

---

☛ show converged plot

Zero-point energy is 3.5905e-01 [eV]

---

2.3.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="400">
               <point coord="1.0     0.0     0.0" label="Gamma"/>
               <point coord="0.625   0.375   0.0" label="K"/>
               <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.0     0.0" label="L"/> 
            </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-band-structure.py (see The python script "PLOT-band-structure.py"), can be used to produce a plot of the phonon-dispersion relations.

$ PLOT-band-structure.py -p -e 0 1400

The resulting PNG file PLOT.png looks like the following.

Converged results as mentioned above can be seen below.

---

☛ show converged plots

Using the path Γ-(K)-X-Γ-L this result is obtained.

A complete visualization of the phonon dispersion of diamond can be achieved using the standard path in reciprocal space for fcc crystals introduced by Setyawan and Curtarolo, Comp. Mat. Sci. 49, 299-312 (2010), which is visualized in the following figure:

The coordinates of the high-symmetry points of the first Brillouin zone belonging to this path are given below.

...     
      <phonondispplot>
         <plot1d>
            <path steps="400">
               <point coord="0.000000 0.000000 0.000000" label="Gamma"/>
               <point coord="0.500000 0.000000 0.500000" label="X"/>
               <point coord="0.500000 0.250000 0.750000" label="W"/>
               <point coord="0.375000 0.375000 0.750000" label="K"/>
               <point coord="0.000000 0.000000 0.000000" label="Gamma"/>
               <point coord="0.500000 0.500000 0.500000" label="L"/>
               <point coord="0.625000 0.250000 0.625000" label="U"/>
               <point coord="0.500000 0.250000 0.750000" label="W"/>
               <point coord="0.500000 0.500000 0.500000" label="L"/>
               <point coord="0.375000 0.375000 0.750000" label="K"/>    
            </path>
         </plot1d>
      </phonondispplot>
...

Using these settings, the following plot is obtained.

---

2.3.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 3$\times$3 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. Phonon properties for cubic boron nitride (cBN)

In this part of the tutorial we will calculate phonons for cubic boron nitride in the zincblende structure to demonstrate the treatment of polar materials. In order to start, create a new working directory, e.g., /home/tutorials/cBN-phonons.

$ cd /home/tutorials
$ mkdir cBN-phonons
$ cd cBN-phonons

In polar materials, a displacement of the atoms creates dipoles which in turn induce a long-ranged electric field. The dipole moment is proportional to an effective charge $Z^\ast_\kappa$ of each atom which also depends on the distribution of the electronic charge and in general is different from the ionic charge $Z_\kappa$. The screening of the electric field is described by the high-frequency dielectric constant $\epsilon^\infty$. In anisotropic materials, both the effective charge and the dielectric constant are tensorial quantities, $Z^\ast_{\alpha\beta,\kappa}$ and $\epsilon^\infty_{\alpha\beta}$. The long-range dipole-fields lead to an additional contribution to the dynamical matrices which is necessary to correctly describe the splitting of longitudinal and transversal optical modes (LO-TO splitting).

---

3.1 Set-up the input file and run calculation

Create the input.xml file with the following content.

<input>
 
   <title>cBN: Phonons</title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="6.850832">
         <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="B.xml" rmt="1.25">
         <atom coord="0.00 0.00 0.00" />
      </species>
 
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.25 0.25 0.25" />
      </species>
 
   </structure>
 
   <groundstate 
      do="fromscratch"
      ngridk="4 4 4"
      swidth="0.0001"
      rgkmax="4.0"
      gmaxvr="14"
      tforce="true">
      <libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
   </groundstate>
 
   <phonons
      do="fromscratch"
      method="dfpt"
      ngridq="2 2 2"
      polar="true">
      <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>

The most important change in comparison to the example of diamond is the attribute polar which we have to set to "true" in order to trigger the calculation of the effective charges $Z^\ast_\kappa$ and dielectric constant $\epsilon^\infty$. The latter is computed using DFPT with an external electric field as perturbation parameter.

Run the calculation using

$ SETUP-excitingroot.sh
$ time exciting_smp &

---

3.2 Output files of a phonon calculation (polar materials)

Additionally to the files discussed in the example of diamond, in polar materials also the files ZSTAR.OUT and EPSINF.OUT are generated, containing the effective charges $Z^\ast_\kappa$ and dielectric constant $\epsilon^\infty$, respectively.

ZSTAR.OUT in this case looks like

# Born effective charge tensors for all atoms.
# Rows correspond to E-field direction.
# Columns correspond to atom displacement direction.
# Acoustic sum rule has been imposed.
#
# species  1 atom   1 (B  1) :      0.000000     0.000000     0.000000
        2.0242834835        0.0000000000        0.0000000000
        0.0000000000        2.0242834835        0.0000000000
        0.0000000000        0.0000000000        2.0242834835
# species  2 atom   1 (N  1) :      0.250000     0.250000     0.250000
       -2.0242834835        0.0000000000        0.0000000000
        0.0000000000       -2.0242834835        0.0000000000
        0.0000000000        0.0000000000       -2.0242834835
# Acoustic sum rule correction (add to each tensor above to get original value)
        0.0005516235        0.0000000000        0.0000000000
        0.0000000000        0.0005516235        0.0000000000
        0.0000000000        0.0000000000        0.0005516235

The acoustic sum rule states that the sum of the charge tensors over all atoms must be zero. By default, we automatically impose this sum rule on the result. The last block displays the correction that was added to the numerical results to fulfill the sum rule.

EPSINF.OUT looks like

# High frequency dielectric tensor (clamped nuclei).
#
        4.5252791835        0.0000000000        0.0000000000
        0.0000000000        4.5252791835        0.0000000000
        0.0000000000        0.0000000000        4.5252791835

In the case of cubic systems such as cBN, all tensors are diagonal and isotropic.

Use again the script CONVERT-au2invcm.py to obtain frequencies at the q-points specified in qpointset.

$ CONVERT-au2invcm.py

############################################

--------------------------------------------
Total number of atoms    :     2
Total number of q-points :     3
--------------------------------------------

q-point   1  is  0.0 0.0 0.0
   mode   1      frequency:     0.0000  cm-1
   mode   2      frequency:     0.0000  cm-1
   mode   3      frequency:     0.0000  cm-1
   mode   4      frequency:   953.6919  cm-1
   mode   5      frequency:   953.6919  cm-1
   mode   6      frequency:   953.6919  cm-1

q-point   2  is  0.5 0.5 0.0
   mode   1      frequency:   707.8988  cm-1
   mode   2      frequency:   707.8988  cm-1
   mode   3      frequency:   826.8165  cm-1
   mode   4      frequency:   826.8165  cm-1
   mode   5      frequency:  1035.6715  cm-1
   mode   6      frequency:  1136.7374  cm-1

q-point   3  is  0.5 0.5 0.5
   mode   1      frequency:   492.4562  cm-1
   mode   2      frequency:   492.4562  cm-1
   mode   3      frequency:   905.2810  cm-1
   mode   4      frequency:   905.2810  cm-1
   mode   5      frequency:   953.5095  cm-1
   mode   6      frequency:  1140.1547  cm-1

############################################

We still see three degenerate optical modes at the Γ-point which is not the correct result. The reason for this is that the polar contribution to the dynamical matrix is not defined for q = 0. Instead, we have to consider the limit lim_q→0 which in general depends on the direction from which we approach the Γ-point. To account for this, choose a point that is slightly displaced from Γ

<input>
 
   <title>cBN: Phonons</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      ...
   </structure>
 
   <groundstate 
      do="skip"
      .../>
   </groundstate>
 
   <phonons
      do="skip"
      ...>
      <qpointset>
         <qpoint> 1d-6 0.0 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.5 </qpoint>
      </qpointset>
   </phonons>
 
</input>

and rerun the calculation (don't forget to set do to "skip" for both groundstate and phonons). You now obtain

$ CONVERT-au2invcm.py

############################################

--------------------------------------------
Total number of atoms    :     2
Total number of q-points :     3
--------------------------------------------

q-point   1  is  1e-06 0.0 0.0
   mode   1      frequency:     0.0005  cm-1
   mode   2      frequency:     0.0005  cm-1
   mode   3      frequency:     2.2034  cm-1
   mode   4      frequency:   953.6919  cm-1
   mode   5      frequency:   953.6919  cm-1
   mode   6      frequency:  1233.9312  cm-1

...

############################################

and we see the splitting between the two transversal and the one longitudinal optical mode. Note also the small frequencies of the acoustic modes slightly away from Γ and the splitting into two transversal and one longitudinal acoustic mode. In cBN, the LO-TO splitting is independent of the direction because the system is isotropic.

---

3.3 Properties from the phonon calculation

Try to repeat the calculation of the phonon density of states and dispersion for cBN. It works exactly the same way as for diamond. You should obtain the following results:

Phonon DOS

Phonon dispersion

---

☛ show converged plot

Phonon DOS

Phonon dispersion

---

4. Useful tips for real-life calculations (DFPT)

4.1 Resume from a previous run

Similar to a ground-state calculation, also a DFPT phonon calculation can be resumed from a previous run. This might prove helpful in cases where an expensive calculation did not finish, e.g., due to given run-time limits or following an unintended crash. In order to restart a calculation from a previous run, set the attribute do to "fromfile" in the phonons element. This will read the density and potential response from the files PHONON_DRHO_fileext_#.OUT and PHONON_DVEFF_fileext_#.OUT and restart the iterative solution of the Sternheimer equation. For that to work, however, it is important that the respective files have been written in the first place. To achieve this, we recommend to make use of the attribute nwrite from the groundstate element in the case of more elaborate calculations. This attribute specifies after how many iterations the current density and potential response are written to file.

Note that subdirectories in which the DYN_fileext_#.OUT files are present will be skipped in following runs.

---

4.2 Parallel execution using multiple processes and MPI

---

☛ Parallelization explained

In general, DFPT phonon calculations can be highly parallelized. There are two levels of parallelization.

1. The outer level parallelizes over pairs (q,I) of q-points and irreps. Each such pair corresponds to one of the subdirectories and constitutes an independent part of the full calculation. The (q,I) pairs do not depend on one another and can be calculated in parallel.
2. The inner level parallelizes over k-points and irrep members d within each of the (q,I) pairs. MPI parallelization over this inner level is also implemented. The product of the number of k-points and irrep members determines the computational load of each of the (q,I) pairs. The computational load is directly proportional to the time needed for a single SCF iteration. It also is the upper limit for the number of processes that can be utilized.

Important information about the computational cost of the full calculation and the parallelization settings that are in use can be found in the file PHONON_RUN_INFO.OUT. Let's discuss the content of the file for the example of diamond above.

We already know the first part from the PHONON_INFO_fileext.OUT file.

================================================================================
= DFPT PHONON CALCULATION INFORMATION                                          =
================================================================================

Irreducible representations will be used.
k-points will be reduced by symmetries.
Response quantities will be symmetrized.

********************************************************************************
* q-points and displacements                                                   *
********************************************************************************

number of q-points:               3
number of phonon modes:           6

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    8        2 (3,3)
   2 ( 0.500000, 0.000000, 0.000000)     12   13        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16   13        3 (2,2,2)

It shows some information about the usage of symmetries and prints all the q-points and irreps with their size (number of irrep members) and the number of (irreducible) k-points for each q-point.

The next part displays information about the computational load of the calculation.

********************************************************************************
* independent parts and computational load                                     *
********************************************************************************

total (remaining) number of independent parts:      9 (9)
total (remaining) computational load:             204 (204)

remaining parts
  ip   iq irrep dim load procs
------------------------------
   1    1     1   3   24     1
   2    1     2   3   24     1
   3    2     1   2   26     1
   4    2     2   1   13     1
   5    2     3   1   13     1
   6    2     4   2   26     1
   7    3     1   2   26     1
   8    3     2   2   26     1
   9    3     3   2   26     1

done parts
  iq irrep dim load
-------------------

In this example, there are 9 independent parts (= number of (q,I) pairs = number of subdirectories). This number equals the sum of the column N_irrep in the table above. The total computational load of the full calculation is 204. I.e., 204 is the largest number of processes that can be effectively utilized in this example.

The two tables below list all the remaining and already completed (none in this example) parts. For each part ip, it gives the index of the q-point iq and the index of the irreducible representation irrep and its dimension dim (= number of the members of irrep). The last two columns give the computational load (= product of number of k-points and irrep members) and the number of processes that will work on this part.

The last part of the file gives further information about the parallelization settings.

********************************************************************************
* parallelization settings                                                     *
********************************************************************************

number of processes:                        1
optimal average load per process:       204.0
minimum load per process:                 204
maximum load per process:                 204
number of idle processes:                   0

proc load (imbalance)
----------------------
   0  204 (+0.0)

In this example, there will be only one process. Both the optimal, minimum and maximum computational load per process is 204'[ in this example and there are no idle processes. The table below shows the load assigned to each process. In parenthesis, the imbalance, i.e., the difference to the optimal average load is provided.

---

Note that one unit of computational load is different for each calculation and it is mainly determined by the number of (L)APW+LO basis functions (= matrix size). Keep reading for a more meaningful example.

---

☛ Choosing parallelization settings

In order to find out how the calculation will be parallelized for a given number of processes you have the option of performing a dry run. Let's assume, we plan to execute the example above with 16 processes. Before we start the actual calculation, we want to find out how the calculation will be parallelized. To this extend, we change the attribute do to "dry" in the phonons element and add the new attribute drynumprocs = "16".

...
   <phonons 
      do="dry" 
      method="dfpt"
      ngridq="2 2 2"
      polar="false"
      drynumprocs="16">
      ...
   </phonons>
...

Running exciting, this case, will only produce the file PHONON_RUN_INFO.OUT. (You might have to delete the DYN_fileext_#.OUT files in the subdiretories if you have already completed the phonon calculation above in the same directory.) This time, the file looks like the following.

================================================================================
= DFPT PHONON CALCULATION INFORMATION                                          =
================================================================================

Irreducible representations will be used.
k-points will be reduced by symmetries.
Response quantities will be symmetrized.

********************************************************************************
* q-points and displacements                                                   *
********************************************************************************

number of q-points:               3
number of phonon modes:           6

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    8        2 (3,3)
   2 ( 0.500000, 0.000000, 0.000000)     12   13        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16   13        3 (2,2,2)

********************************************************************************
* independent parts and computational load                                     *
********************************************************************************

total (remaining) number of independent parts:      9 (9)
total (remaining) computational load:             204 (204)

remaining parts
  ip   iq irrep dim load procs
------------------------------
   1    1     1   3   24     3
   2    1     2   3   24    12
   3    2     1   2   26    13
   4    2     2   1   13    13
   5    2     3   1   13    13
   6    2     4   2   26    13
   7    3     1   2   26    13
   8    3     2   2   26    13
   9    3     3   2   26    13

done parts
  iq irrep dim load
-------------------

********************************************************************************
* parallelization settings                                                     *
********************************************************************************

number of processes:                       16
optimal average load per process:        12.8
minimum load per process:                   8
maximum load per process:                  14
number of idle processes:                   0

proc load (imbalance)
----------------------
   0   14 (+1.2)
   1   14 (+1.2)
   2   14 (+1.2)
   3   14 (+1.2)
   4   14 (+1.2)
   5   14 (+1.2)
   6   14 (+1.2)
   7   14 (+1.2)
   8   14 (+1.2)
   9   14 (+1.2)
  10   14 (+1.2)
  11   14 (+1.2)
  12   12 (-0.8)
  13    8 (-4.8)
  14    8 (-4.8)
  15    8 (-4.8)

The table in the second section of the file now shows how many processes will be working on which part (e.g., 12 processes will work on part 2). From the third section we can tell that 16 is a reasonable number of processes for this example, since the optimal average load is 12.8 and the maximum load per process is 14 which is the one more than the smallest possible integer greater than the average load and meaning only little "waste" of resources. You will have bad parallelization settings if you either have a large number of idle processes or a very uneven distribution of load among the processes, i.e., the maximum load per process being much larger than the optimal average load.

---

5. Alternative: Phonons from supercell calculations

Alternatively to using DFPT to compute phonons, you can also calculate phonons from supercell calculations. In that case, for each q-point a corresponding supercell is constructed and each equivalent atom is displaced along the three Cartesian directions by a small amount. The actual displacement pattern of these atoms is a cosinus or sinus wave with periodicity defined by the wavevector q. The amplitude of the wavelike displacement pattern is defined by the value Δ_ph of the attribute deltaph. Finally, elements of the dynamical matrix are obtained using a finite difference procedure from atomic forces calculated for displacement patterns with amplitude Δ_ph and 2Δ_ph (-Δ_ph and Δ_ph at q=0). In general, using DFPT leads to much faster calculations than using the supercell approach. However, supercell phonons can still be very useful to investigate anharmonic effects beyond the harmonic approximation (anharmonicity). Further, not all features available in ground-state calculations are implemented within DFPT, but they can be used using supercells. Examples are van-der-Waals interactions or spin-orbit coupling.

---

5.1 Set up the supercell phonon calculation

In order to trigger a supercell calculation, we have switch the attribute method to "sc" and set the displacement amplitude (in Bohr units) using deltaph.

...
   <phonons 
      do="fromscratch" 
      method="sc"
      ngridq="2 2 2"
      deltaph="0.03"
      polar="false">
      ...
   </phonons>
...

---

5.2 Output files of a phonon calculation

The name of the subdirectories created from a supercell calculation is different from the DFPT case. The names of the directories include information, in order of appearance, on the q vector (Q…), the species of the displaced atoms (S…), the index in the unperturbed primitive cell of the displaced atoms (A…), and the polarization direction of the displacement (P…). Example of such names are Q0000_0000_0000_S01_A001_P1, Q0102_0000_0000_S01_A001_P1, etc.

Inside the subdirectories you will find the files INFO_fileext.OUT which are standard INFO.OUT files of the various groundstate calculations performed in the finite difference scheme. They replace the PHONON_INFO_fileext.OUT files from the DFPT case.

---

5.3 Polar materials from supercell calculations

A supercell calculation of Γ-point phonons also calculates the Born effective charge tensors $Z^\ast_\kappa$ and generates the file ZSTAR.OUT. The dielectric constant $\epsilon^\infty$ in the corresponding file EPSINF.OUT, however, is not computed and must be provided manually. Once both files are present, all further phonon properties will correctly account for the polar character if you set polar to "true".

---

5.4 Phonon properties for diamond

As an exercise try to repeat the phonon calculation for diamond but now using the supercell method. Here, we will only present the results you should obtain using a 4$\times$4$\times$4 k-grid, rgkmax = "4.0", and 2 q-points in each direction.

Phonon frequencies

$ CONVERT-au2invcm.py

############################################

--------------------------------------------
Total number of atoms    :     2
Total number of q-points :     3
--------------------------------------------

q-point   1  is  0.0 0.0 0.0
   mode   1      frequency:     0.0000  cm-1
   mode   2      frequency:     0.0000  cm-1
   mode   3      frequency:     0.0000  cm-1
   mode   4      frequency:  1258.2158  cm-1
   mode   5      frequency:  1258.2158  cm-1
   mode   6      frequency:  1258.2158  cm-1

q-point   2  is  0.5 0.5 0.0
   mode   1      frequency:   777.4162  cm-1
   mode   2      frequency:   777.4162  cm-1
   mode   3      frequency:  1014.2525  cm-1
   mode   4      frequency:  1014.2525  cm-1
   mode   5      frequency:  1195.8142  cm-1
   mode   6      frequency:  1195.8142  cm-1

q-point   3  is  0.5 0.5 0.5
   mode   1      frequency:   540.8468  cm-1
   mode   2      frequency:   540.8468  cm-1
   mode   3      frequency:  1009.7927  cm-1
   mode   4      frequency:  1153.6611  cm-1
   mode   5      frequency:  1153.6611  cm-1
   mode   6      frequency:  1267.8193  cm-1

############################################

Note the differences in comparison with the DFPT result. These differences become much smaller for more converged numerical settings.

Phonon DOS

Phonon dispersion