by Pasquale Pavone for exciting carbon
Purpose: In this tutorial, you will learn how to set up and execute a series of calculations for a crystal in the diamond structure, where the atoms are displaced according to a phonon pattern with the periodicity of the X point in the Brillouin zone. Additionally, it will be explained how to obtain the phonon frequency of TA, LA, LO, and TO modes at X, by calculating the derivatives of the energyvsdisplacement curves at zero displacement.
Table of Contents

0. Define relevant environment variables
Proceed in the same way as for Phonons at Γ in diamondstructure crystals.
Here, is a list of the scripts which are relevant for this tutorial with a short description.
 SETUPdiamondphononx.py: Python script for generating structures with displaced atoms.
 EXECUTEdiamondphononx.sh: (Bash) shell script for running a series of exciting calculations.
 CHECKFITenergyvsdisplacement.py: Python script for extracting the derivatives at zero displacement of the energyvsdisplacement curves.
 PLOTenergy.py: Python visualization tool for the energyvsdisplacement curves.
 PLOTstatus.py: Python visualization tool for the RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
 PLOTcheckderiv.py: Python visualization tool for the fit of the derivatives of the energyvsdisplacement curves at zero displacement.
1. Set up the calculations
i) Preparation of the input file
The first step is to create a directory for each system that you want to investigate. In this tutorial, we consider as an example the calculation of the energyvsdisplacement curves for carbon in the diamond structure. Thus, we will create a directory diamondphononx and we move inside it.
$ mkdir diamondphononx
$ cd diamondphononx
Inside this directory, we copy the file input.xml displayed in Phonons at Γ in diamondstructure crystals. Please, remember that the input file for an exciting calculation must always be called input.xml. After creating the input file, set the path for the species file using the command
$ SETUPexcitingroot.sh
ii) Periodicity of phonons at the X point
The perturbed crystal, formed by displacing the atoms according to a phonon pattern at the X point in the Brillouin zone, has a unit cell which size is doubled with respect to the unperturbed crystal. The new unit cell has four atoms and a tetragonal Bravais lattice with basis vectors
 (1, 0, 0) a / √2;
 (0, 1, 0) a / √2;
 (0, 0, 1) a;
where a is the cubic lattice constant. Furthermore, the positions of the four atoms in the doubled unit cell are (notice that for tetragonal structures all crystal axes are parallel to Cartesian directions):
 (0, 0, 0);
 (0, 1/2, 1/4);
 (1/2, 1/2, 1/2);
 (1/2, 1, 3/4).
According to this, for performing the perturbed calculations, the input file showed above should be changed to take into account the new structure. This is done by the script SETUPdiamondphononx.py described in the next subsection. The script will generate input files with the new geometry. For instance, for vanishing displacement, the crystal and species elements of the new input look like this:
... <crystal scale=" 4.7707080313"> <basevect> 1.0000000000 0.0000000000 0.0000000000 </basevect> <basevect> 0.0000000000 1.0000000000 0.0000000000 </basevect> <basevect> 0.0000000000 0.0000000000 1.4142135624 </basevect> </crystal> <species speciesfile="C.xml" rmt="1.25"> <atom coord=" 0.0000000000 0.0000000000 0.0000000000"/> <atom coord=" 0.0000000000 0.5000000000 0.2500000000"/> <atom coord=" 0.5000000000 0.5000000000 0.5000000000"/> <atom coord=" 0.5000000000 1.0000000000 0.7500000000"/> </species> ...
Furthermore, when setting up a new calculation for a tetragonal system, the values defined by the attribute ngridk in the element groundstate will be also changed according to the new geometry. The new values are chosen in such a way to keep the same density of kpoints for both the unperturbed unit cell and the tetragonal supercell.
Note that you do NOT need to modify the file input.xml by adding explicitly the configuration with 4 atoms in the super cell. The script SETUPdiamondphononx.py described in the next section will do it for you.
iii) Generation of input files for the different structures
In order to generate input files for a series of different structures you have to run the script SETUPdiamondphononx.py. Notice that the script SETUPdiamondphononx.py always generates a working directory containing input files for different atomic displacements. Results of the current calculations will be also stored in the working directory. The directory name can be specified by adding the name in the command line.
$ SETUPdiamondphononx.py DIRECTORYNAME
If no name is given, the script use the default name workdir. Very important: The working directory is overwritten each time you execute the script SETUPdiamondphononx.py. Therefore, choose different names for different calculations.
The script SETUPdiamondphononx.py produces the following output on the screen (using xphononta as working directory).
$ SETUPdiamondphononx.py xphononta
Enter maximum displacement umax [u/(lat. parameter)] >>>> 0.035
Enter the number of displacements in [umax,umax] >>>> 21

List of labels for phonon modes at the X point in the diamond structure

TA => Transverse acoustic mode
LA => Longitudinal acoustic mode
TO => Transverse optical mode
LO => Longitudinal optical mode

Enter label for phonon mode [1 choice] >>>> TA
$
In this example, the (on screen) input entries are preceded by the symbol ">>>>". The entry values must be typed on the screen when requested. The first entry (in our example 0.035) represents the absolute value of the maximum displacement (for each component, in crystal coordinates) for which we want to perform the calculation. The second entry (21) is the number of structures equally spaced in the displacement of the second atom, which are generated between the maximum negative displacement and the maximum positive one. The third (last) entry (TA) is the code of the phonon mode that one wants to investigate.
2. Execute the calculations
To execute the series of calculation with input files created by SETUPdiamondphononx.py you have to run the script EXECUTEdiamondphononx.sh. If a name for the working directory has been specified, then you must give it here, too.
$ EXECUTEdiamondphononx.sh xphononta
Running exciting for file input11.xml 
...
Run completed for file input21.xml 
$
Notice that, due to the symmetry of the phonon pattern at X, the script executes only calculations with input files corresponding to non negative displacements. This is the reason why the first executed calculation is for the file input11.xml instead of for input01.xml.
After the complete run, inside the directory xphononta the results of the calculation for the input file inputi.xml are contained in the subdirectory rundiri where i is running from 11 to the total number of strain values. Data for the energyvsdisplacement curves are contained in the files energyvsdisplacement.
3. Postprocessing: Extract energy derivatives
In order to analyse the results of the calculations, you first have to move to the working directory.
$ cd xphononta
At this point, you can use the python script CHECKFITenergyvsdisplacement.py for extracting the derivatives of the energyvsdisplacement curves at zero displacement.
$ CHECKFITenergyvsdisplacement.py
Enter maximum displacement for the fit >>>> 0.035
Enter the order of derivative >>>> 2
Enter atomic mass in [amu] >>>> 12.01
_____________________________________________
Xphononcalculation
#############################################
Fit data
Maximum value of the displacement ==> 0.035
Number of displacement values used ==> 21
Fit results for the derivative of order 2
Polynomial of order 2 ==> 767.82 [cm1]
Polynomial of order 3 ==> 767.82 [cm1]
Polynomial of order 4 ==> 738.11 [cm1]
Polynomial of order 5 ==> 738.11 [cm1]
Polynomial of order 6 ==> 729.95 [cm1]
Polynomial of order 7 ==> 729.95 [cm1]
Polynomial of order 2 ==> 23.0188 [THz]
Polynomial of order 3 ==> 23.0188 [THz]
Polynomial of order 4 ==> 22.1280 [THz]
Polynomial of order 5 ==> 22.1280 [THz]
Polynomial of order 6 ==> 21.8835 [THz]
Polynomial of order 7 ==> 21.8835 [THz]
#############################################
$
In this example, the input entries are preceded by the symbol ">>>>". The entry values must be typed on the screen when requested. The first entry (in our example 0.035) represents the absolute value of the maximum displacement for which we want to perform the analysis. The second entry (2) is the order of the derivative that we want to obtain. Finally, the third entry (12.01) is the atomic mass of carbon in units of [amu].
NOTICE that, in this example, the values which are given above as output on the screen are not directly the actual second derivative of the energyvsdisplacement curves, but the values of the frequency, i.e., combinations involving the square root of the derivative (see Section 5 for the explanation).
The script CHECKFITenergyvsdisplacement.py generates the output files checkenergyderivatives and orderofderivative, which can be used in the postprocessing analysis.
4. Postprocessing: Visualization tools
All scripts, except PLOTatomforce.py, which are used as visualization tools of in Phonons at Γ in diamondstructure crystals, can be used for this tutorial, too. In the following, you can find some example of plot obtained with these tools.
PLOTenergy.py
PLOTcheckderiv.py
PLOTstatus.py rundir21
5. Postprocessing: How to derive the optical phonon frequency at the X point
The total energy per unit cell of a crystal where the atoms are displaced by their equilibrium positions can be written as a Taylor series
E({u}) = E_{0} + 1/2 Σ_{R,R'} u(R) Φ(R,R') u(R') + O(u^{3})
where E_{0} is the energy corresponding to the equilibrium structure, u(R) is an atomic displacement in the cell R, Φ(R,R') is the interatomic forceconstants matrix, and O(u^{3}) includes all orders higher than two in the displacements u. For small displacements, terms beyond the harmonic approximation (which retains only terms that are quadratic in u) can be neglected. For the diamond structure the displacements of the 4 atoms in the "doubled" unit cell corresponding to the phonon modes at the X point can be given as (a is the cubic lattice constant):
TA phonon
 u_{1} = a (0, u, 0)
 u_{2} = a (0, u, 0)
 u_{3} = a (0, u, 0)
 u_{4} = a (0, u, 0)
LA phonon
 u_{1} = a (0, 0, u)
 u_{2} = a (0, 0, u)
 u_{3} = a (0, 0, u)
 u_{4} = a (0, 0, u)
TO phonon
 u_{1} = a (0, u, 0)
 u_{2} = a (0, u, 0)
 u_{3} = a (0, u, 0)
 u_{4} = a (0, u, 0)
LO phonon
 u_{1} = a (0, 0, u)
 u_{2} = a (0, 0, u)
 u_{3} = a (0, 0, u)
 u_{4} = a (0, 0, u)
Notice, that in the diamond structure the LA and LO phonons are degenerate (they have the same frequency).
For a phonon mode at the X point, the total energy per unit cell (containing 2 atoms) of the crystal can be written as (see details here):
E({u}) = E_{0} + m a^{2} ω(X)^{2} u^{2}
Using appropriate values for m and a and their units, the phonon frequency is given directly from the square root of the second derivative of the energyvsdisplacement curves. This is the output of the fit procedure exposed above if the order of the derivative is taken equal to 2. The optical phonon frequency is obtained from the plateau value of the secondorder derivative at zero displacement of the energyvsdisplacement curves. For more details on the procedure used for extracting numerical derivatives, see Energy vs. strain calculations, where the same approach is used for calculating elastic constants.
Exercises
 The frequency values which are obtained in this example using the quoted parameters, are not realistic. Check the convergence of the phonon frequencies at the X point for diamond with respect of the value of the ngridk and rgkmax attributes using derivatives of the total energy.
 Calculate all modes at the X point
 TA mode
 TO mode
 LA mode
 LO mode
 Calculate all frequency modes for different values of the lattice constant.
 Repeat all the calculation for silicon.