by Sebastian Tillack for exciting fluorine
Purpose: In this tutorial, you will learn how to perform a phonon calculation using densityfunctional perturbation theory (DFPT). You will work through the example of diamond.
Note: We strongly recommend to do the supercell phonon tutorial Phonon properties of diamondstructure crystals (SuperCell) first to familiarize yourself with the standard phonon input parameters. All phononderived properties from the supercell tutorial can be applied in the same way to DFPT phonons.
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. We further assume that you have already completed Phonon properties of diamondstructure crystals (SuperCell).
In order to start, create a new working directory, e.g., /home/tutorials/diamondphononsdfpt.
$ cd /home/tutorials
$ mkdir diamondphononsdfpt
$ cd diamondphononsdfpt
1. Theoretical Background
2. Running a DFPT phonon calculation
In comparison to the input file from Phonon properties of diamondstructure crystals (SuperCell) two essential changes in the input file are necessary.
First, in the groundstate element, 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. In the case of the PBE xcfunctional, the corresponding line in the input file reads
<libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
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.
The second change is to switch the method used to compute phonons to "dfpt" to trigger a DFPT calculation.
... <phonons do="fromscratch" method="dfpt" ngridq="2 2 2"> ... </phonons> ...
The complete input file for this calculation is given in the following.
<input> <title>Diamond: DFPT Phonons</title> <structure speciespath="$EXCITINGROOT/species"> <crystal scale="6.7468"> <basevect>0.5 0.5 0.0</basevect> <basevect>0.5 0.0 0.5</basevect> <basevect>0.0 0.5 0.5</basevect> </crystal> <species speciesfile="C.xml" rmt="1.25"> <atom coord="0.00 0.00 0.00" /> <atom coord="0.25 0.25 0.25" /> </species> </structure> <groundstate do="fromscratch" ngridk="2 2 2" swidth="0.0001" rgkmax="4.0" gmaxvr="14" maxscl="25" > <libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/> </groundstate> <phonons do="fromscratch" method="dfpt" ngridq="2 2 2"> <qpointset> <qpoint> 0.0 0.0 0.0 </qpoint> <qpoint> 0.5 0.5 0.0 </qpoint> <qpoint> 0.5 0.5 0.5 </qpoint> </qpointset> </phonons> </input>
Run the calculation as usual with
$ SETUPexcitingroot.sh
$ time exciting_smp &
3. Output of a DFPT phonon calculation
3.1 Common output files
The most important output files such as PHONON.OUT and the dynamical matrices in DYN_Q####_####_####_S##_A###_P#.OUT are the same as in Phonon properties of diamondstructure crystals (SuperCell). Hence, you still can see the resulting phonon frequencies using
$ 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: 1537.2055 cm1
mode 5 frequency: 1537.2055 cm1
mode 6 frequency: 1537.2055 cm1
qpoint 2 is 0.5 0.5 0.0
mode 1 frequency: 677.4866 cm1
mode 2 frequency: 677.4866 cm1
mode 3 frequency: 1005.5120 cm1
mode 4 frequency: 1005.5120 cm1
mode 5 frequency: 1202.3322 cm1
mode 6 frequency: 1202.3322 cm1
qpoint 3 is 0.5 0.5 0.5
mode 1 frequency: 452.2903 cm1
mode 2 frequency: 452.2903 cm1
mode 3 frequency: 894.6191 cm1
mode 4 frequency: 1290.8184 cm1
mode 5 frequency: 1290.8184 cm1
mode 6 frequency: 1377.7608 cm1
############################################
Note the differences in the resulting frequencies compared to the supercell approach. These are due to the different methods and become much smaller for more converged numerical settings.
In order to get the phonon density of states and the full phonon dispersion curves, one can proceed as in tutorial Phonon properties of diamondstructure crystals (SuperCell). The resulting phonon densityofstates will look like:
The obtained phonon dispersion curves are shown below.
3.2 DFPT specific output
A DFPT phonon calculation produces different subdirectories than a supercell calculation. In a DFPT calculation we use special symmetry adapted displacement patterns, the irreps, instead of the canonical displacements of each atom in each of the three Cartesian directions as done in the supercell calculation. Hence, the names of the subdirectories in a DFPT calculation are Q####_####_####_I###, where the first part specifies the qpoint in the grid and the second part labels the irrep at that specific qpoint. 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.
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 3 2 (3*,3)
2 ( 0.500000, 0.000000, 0.000000) 12 4 4 (2,1,1,2)
3 ( 0.500000, 0.500000, 0.000000) 16 4 3 (2,2,2)
current task / total (remaining) number of tasks: 1/ 9 (9)
load of current task / total (remaining) load: 9/ 66 (66)
procs for current task / total number of procs: 1/ 1
average / maximum load per proc: 66/ 66
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 4.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 *
********************************************************************************
Initial density response read from file.
Sternheimer SCF loop iteration 1
maximum change in potential response: 0.101255
Fermi energy response: 0.00000
Time spent for iteration (seconds): 0.77
Sternheimer SCF loop iteration 2
maximum change in potential response: 0.697163E01
Fermi energy response: 0.00000
Time spent for iteration (seconds): 0.72
...
4. Useful tips for reallife calculations
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.