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 105 cm-1. ### 1. Theoretical Background ### 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)

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. ###### 2.3.2 Thermodynamic properties The zero-point energy, Ezp, and further thermodynamic properties can be computed from the phonon DOS g$(\omega)$(see Input reference for further details). 1. the vibrational free energy Fvib= Uvib -TSvib, 2. the vibrational internal energy Uvib, 3. the entropic contribution to the vibrational free energy TSvib, 4. the vibrational entropy Svib 5. the heat capacity Cv. 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.

###### 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. ###### 2.3.4 Other useful features ### 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 limq→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

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

### 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 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