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 phonondispersion relations. You will work through the examples of diamond and cubic Boron nitride.
Table of Contents

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.
 FREEfromdos.py: Python script for calculating thermodynamic properties.
 PLOTdos.py: Python visualization tool for the phonon density of states.
 PLOTbandstructure.py: Python visualization tool for phonondispersion curves.
 PLOTfiles.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^{5} 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/diamondphonons.
$ cd /home/tutorials
$ mkdir diamondphonons
$ cd diamondphonons
2.1 Setup 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 qpoints mesh is defined in the same way as the kpoints 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 qgrid 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 exchangecorrelation potential from the libxc library instead of the exciting internal implementation. The reason for this change is that a DFPT calculation requires the xckernel (the second derivatives of the xcenergy) which is not natively implemented in exciting and is always obtained from the libxc library. For full consistency, we recommend to also use the xcpotential (first derivative of the xcenergy) from libxc. However, in principle you can also use the exciting internal implementations of the xcpotential using the attribute xctype. In this case, only the xckernel 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 qpoints (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 qpoints specified in each subelement qpoint. In the example above, calculations will be performed for Γ, X, and L phonons (in the standard notation for highsymmetry points in the Brillouin zone of fcc crystals). Note that qpoints 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
$ SETUPexcitingroot.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 qpoint in the grid and the second part labels the irrep at that specific qpoint. Irreps are special symmetry adapted displacement patterns. Inside of each subdirectory, you will find the following files:

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.
kpoints are reduced by symmetries.
Response quantities are symmetrized.
qpoints and symmetries:
iq qvector (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) qpoints 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 kpoints in the irreducible Brillouin zone obtained with the symmetries in the small group. N_irrep denotes the number of irreps for the given qpoint and the values in parenthesis give their respective size (or equivalently the number of members per irrep). The qpoint 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 selfconsistency 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.960508E01
Fermi energy response: 0.00000
Time spent for iteration (seconds): 1.20
Sternheimer SCF loop iteration 2
maximum change in potential response: 0.646932E01
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 qpoints 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 CONVERTau2invcm.py to obtain frequencies expressed in units of cm^{1}.
$ CONVERTau2invcm.py
############################################

Total number of atoms : 2
Total number of qpoints : 3

qpoint 1 is 0.0 0.0 0.0
mode 1 frequency: 0.0000 cm1
mode 2 frequency: 0.0000 cm1
mode 3 frequency: 0.0000 cm1
mode 4 frequency: 1244.0005 cm1
mode 5 frequency: 1244.0005 cm1
mode 6 frequency: 1244.0005 cm1
qpoint 2 is 0.5 0.5 0.0
mode 1 frequency: 782.9265 cm1
mode 2 frequency: 782.9265 cm1
mode 3 frequency: 1006.6311 cm1
mode 4 frequency: 1006.6311 cm1
mode 5 frequency: 1195.6297 cm1
mode 6 frequency: 1195.6297 cm1
qpoint 3 is 0.5 0.5 0.5
mode 1 frequency: 539.9045 cm1
mode 2 frequency: 539.9045 cm1
mode 3 frequency: 998.7204 cm1
mode 4 frequency: 1149.7671 cm1
mode 5 frequency: 1149.7671 cm1
mode 6 frequency: 1275.9419 cm1
############################################
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
$ SETUPexcitingroot.sh
$ time exciting_smp &
Data in PHDOS.OUT can be plotted directly. You can use the script PLOTdos.py (see The python script "PLOTdos.py"), by typing on the command line
$ PLOTdos.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 kgrid, rgkmax = "7.0", and 4 qpoints in each direction) can be found below.
2.3.2 Thermodynamic properties
The zeropoint energy, E_{zp}, and further thermodynamic properties can be computed from the phonon DOS g$(\omega)$ (see Input reference for further details).
 the vibrational free energy F_{vib}= U_{vib} TS_{vib},
 the vibrational internal energy U_{vib},
 the entropic contribution to the vibrational free energy TS_{vib},
 the vibrational entropy S_{vib}
 the heat capacity C_{v}.
In order to calculate and visualize these thermodynamic properties, you can use the script FREEfromdos.py as follows
$ FREEfromdos.py 0 1400 160 eV
Zeropoint energy is 3.5124e01 [eV]
Created PNG output "PLOTfreeenergies.png"
Created PNG output "PLOTentropy.png"
Created PNG output "PLOTheatcapacity.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 zeropoint energy is written on the screen output. The script also generates three PNG output files: PLOTfreeenergies.png, PLOTentropy.png, and PLOTheatcapacity.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 Phonondispersion relations
Phonondispersion 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 PLOTbandstructure.py (see The python script "PLOTbandstructure.py"), can be used to produce a plot of the phonondispersion relations.
$ PLOTbandstructure.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/cBNphonons.
$ cd /home/tutorials
$ mkdir cBNphonons
$ cd cBNphonons
In polar materials, a displacement of the atoms creates dipoles which in turn induce a longranged 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 highfrequency 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 longrange dipolefields lead to an additional contribution to the dynamical matrices which is necessary to correctly describe the splitting of longitudinal and transversal optical modes (LOTO splitting).
3.1 Setup 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
$ SETUPexcitingroot.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 Efield 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 CONVERTau2invcm.py to obtain frequencies at the qpoints specified in qpointset.
$ CONVERTau2invcm.py
############################################

Total number of atoms : 2
Total number of qpoints : 3

qpoint 1 is 0.0 0.0 0.0
mode 1 frequency: 0.0000 cm1
mode 2 frequency: 0.0000 cm1
mode 3 frequency: 0.0000 cm1
mode 4 frequency: 953.6919 cm1
mode 5 frequency: 953.6919 cm1
mode 6 frequency: 953.6919 cm1
qpoint 2 is 0.5 0.5 0.0
mode 1 frequency: 707.8988 cm1
mode 2 frequency: 707.8988 cm1
mode 3 frequency: 826.8165 cm1
mode 4 frequency: 826.8165 cm1
mode 5 frequency: 1035.6715 cm1
mode 6 frequency: 1136.7374 cm1
qpoint 3 is 0.5 0.5 0.5
mode 1 frequency: 492.4562 cm1
mode 2 frequency: 492.4562 cm1
mode 3 frequency: 905.2810 cm1
mode 4 frequency: 905.2810 cm1
mode 5 frequency: 953.5095 cm1
mode 6 frequency: 1140.1547 cm1
############################################
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> 1d6 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
$ CONVERTau2invcm.py
############################################

Total number of atoms : 2
Total number of qpoints : 3

qpoint 1 is 1e06 0.0 0.0
mode 1 frequency: 0.0005 cm1
mode 2 frequency: 0.0005 cm1
mode 3 frequency: 2.2034 cm1
mode 4 frequency: 953.6919 cm1
mode 5 frequency: 953.6919 cm1
mode 6 frequency: 1233.9312 cm1
...
############################################
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 LOTO 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 reallife calculations (DFPT)
4.1 Resume from a previous run
Similar to a groundstate 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 runtime 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
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 qpoint 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 groundstate calculations are implemented within DFPT, but they can be used using supercells. Examples are vanderWaals interactions or spinorbit 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 kgrid, rgkmax = "4.0", and 2 qpoints in each direction.
Phonon frequencies
$ CONVERTau2invcm.py
############################################

Total number of atoms : 2
Total number of qpoints : 3

qpoint 1 is 0.0 0.0 0.0
mode 1 frequency: 0.0000 cm1
mode 2 frequency: 0.0000 cm1
mode 3 frequency: 0.0000 cm1
mode 4 frequency: 1258.2158 cm1
mode 5 frequency: 1258.2158 cm1
mode 6 frequency: 1258.2158 cm1
qpoint 2 is 0.5 0.5 0.0
mode 1 frequency: 777.4162 cm1
mode 2 frequency: 777.4162 cm1
mode 3 frequency: 1014.2525 cm1
mode 4 frequency: 1014.2525 cm1
mode 5 frequency: 1195.8142 cm1
mode 6 frequency: 1195.8142 cm1
qpoint 3 is 0.5 0.5 0.5
mode 1 frequency: 540.8468 cm1
mode 2 frequency: 540.8468 cm1
mode 3 frequency: 1009.7927 cm1
mode 4 frequency: 1153.6611 cm1
mode 5 frequency: 1153.6611 cm1
mode 6 frequency: 1267.8193 cm1
############################################
Note the differences in comparison with the DFPT result. These differences become much smaller for more converged numerical settings.
Phonon DOS
Phonon dispersion