by Pasquale Pavone for exciting oxygen
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 second atom in the unit cell is displaced along the cube diagonal. Additionally, it will be explained how to obtain the phonon frequency of the optical modes at q=0 (Γ point), by calculating the derivatives of the energyvsdisplacement and forcevsdisplacement curves at zero displacement.
Table of Contents

0. Define relevant environment variables
Read the following paragraphs before starting with the rest of this tutorial!
Before starting, be sure that relevant environment variables are already defined as specified in 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.
 SETUPdiamondphonong.py: Python script for generating structures with displaced atoms.
 EXECUTEdiamondphonong.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.
 CHECKFITforcevsdisplacement.py: Python script for extracting the derivatives at zero displacement of the forcevsdisplacement curves.
 PLOTenergy.py: Python visualization tool for the energyvsdisplacement curves.
 PLOTatomforce.py: Python visualization tool for the forcevsdisplacement 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.
From now on the symbol $ will indicate the shell prompt.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries. Gnuplot (visualization tools).
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 diamondphonong and we move inside it.
$ mkdir diamondphonong
$ cd diamondphonong
Inside this directory, we create (or copy from a previous calculation) the file input.xml corresponding to a calculation for the equilibrium structure of diamond. This file could look like the following
<input> <title>Diamond: Phonons at Gamma</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 ngridk="2 2 2" xctype="GGA_PBE" swidth="0.0001" rgkmax="4.0" gmaxvr="14" tforce="true"> </groundstate> </input>
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
Furthermore
 do not perform any structural optimization.
 perform the calculation of the force acting on the atoms by adding the tforce = "true" attribute to the groundstate element:
... <groundstate ... tforce="true"> </groundstate> ...
ii) 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 SETUPdiamondphonong.py. Notice that the script SETUPdiamondphonong.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 (in this example gamma1) can be specified by adding the name in the command line.
$ SETUPdiamondphonong.py gamma1
If the directory name is not given explicitly, the script SETUPdiamondphonong.py generates a directory called by default workdir.
Very important: The working directory is overwritten each time you execute the script SETUPdiamondphonong.py. Therefore, choose different names for different calculations.
The script SETUPdiamondphonong.py produces the following output on the screen (using gamma1 as working directory).
$ SETUPdiamondphonong.py gamma1
Enter maximum displacement umax [u/(lat. parameter)] >>>> 0.025
Enter the number of displacements in [umax,umax] >>>> 11
$
In this example, (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.025) 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 (11) 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.
2. Execute the calculations
To execute the series of calculations with input files created by SETUPdiamondphonong.py you have to run the script EXECUTEdiamondphonong.sh. If a name for the working directory has been specified, then you must give it here, too.
$ EXECUTEdiamondphonong.sh gamma1
Running exciting for file input01.xml 
...
Run completed for file input11.xml 
$
After the complete run, inside the directory gamma1 the results of the calculation for the input file inputi.xml are contained in the subdirectory rundiri where i is running from 01 to the total number of strain values. The data for the energyvsdisplacement and forcevsdisplacement curves are contained in the files energyvsdisplacement and forcevsdisplacement, respectively. More precisely, the file forcevsdisplacement contains the x component of the force acting on the second atom in the unit cell as a function of the displacement.
3. Postprocessing: Extract energy derivatives
In order to analyse the results of the calculations, you first have to move to the working directory.
$ cd gamma1
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.025
Enter the order of derivative >>>> 2
Enter atomic mass in [amu] >>>> 12.01
#############################################
Fit data
Maximum value of the displacement ==> 0.025
Number of displacement values used ==> 11
Fit results for the derivative of order 2
Polynomial of order 2 ==> 1602.12 [cm1]
Polynomial of order 3 ==> 1602.12 [cm1]
Polynomial of order 4 ==> 1584.17 [cm1]
Polynomial of order 5 ==> 1584.17 [cm1]
Polynomial of order 6 ==> 1584.05 [cm1]
Polynomial of order 7 ==> 1584.05 [cm1]
Polynomial of order 2 ==> 48.0303 [THz]
Polynomial of order 3 ==> 48.0303 [THz]
Polynomial of order 4 ==> 47.4923 [THz]
Polynomial of order 5 ==> 47.4923 [THz]
Polynomial of order 6 ==> 47.4885 [THz]
Polynomial of order 7 ==> 47.4885 [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.025) 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 6. for the explanation).
The script CHECKFITenergyvsdisplacement.py generates the output files checkenergyderivatives and orderofderivative, which can be used in the postprocessing analysis.
4. Postprocessing: Extract force derivatives
The procedure for extracting the derivatives of the forcevsdisplacement curves is very similar to the one exposed in the previous Section for the energyvsdisplacement curves. You can use the python script CHECKFITforcevsdisplacement.py for extracting the derivatives of the forcevsdisplacement curves at zero displacement. The meaning of the input entries for this script is the same as for CHECKFITenergyvsdisplacement.py (see above). The only difference is that the derivative which one has to calculate in order to extract the phonon frequencies is the first one. Here, is an example of execution of the script CHECKFITforcevsdisplacement.py.
$ CHECKFITforcevsdisplacement.py
Enter maximum displacement for the fit >>>> 0.025
Enter the order of derivative >>>> 1
Enter atomic mass >>>> 12.01
#############################################
Fit data
Maximum value of the displacement ==> 0.025
Number of displacement values used ==> 11
Fit results for the derivative of order 1
Polynomial of order 1 ==> 1605.97 [cm1]
Polynomial of order 2 ==> 1605.97 [cm1]
Polynomial of order 3 ==> 1580.60 [cm1]
Polynomial of order 4 ==> 1580.60 [cm1]
Polynomial of order 5 ==> 1580.62 [cm1]
Polynomial of order 6 ==> 1580.62 [cm1]
Polynomial of order 1 ==> 48.1456 [THz]
Polynomial of order 2 ==> 48.1456 [THz]
Polynomial of order 3 ==> 47.3852 [THz]
Polynomial of order 4 ==> 47.3852 [THz]
Polynomial of order 5 ==> 47.3859 [THz]
Polynomial of order 6 ==> 47.3859 [THz]
#############################################
$
5. Postprocessing: Visualization tools
All the scripts mentioned here must be executed in the directory where the energyvsdisplacement, forcevsdisplacement, checkenergyderivatives, and orderofderivative files are located. The scripts produce as output a PostScript file named PLOT.ps .
i) PLOTenergy.py
This script allows for the visualization of the energyvsdisplacement curve. It is executed as follows.
$ PLOTenergy.py
ii) PLOTcheckderiv.py
This is a very important tool that allows to represent the dependence of the calculated derivatives of the energyvsdisplacement and forcevsdisplacement curves on
 the range of points included in the fitting procedure ("maximum displacement u"),
 the maximum degree of the polynomial used in the fitting procedure ("n").
The script PLOTcheckderiv.py requires as input the checkenergyderivatives and orderofderivative files generated by CHECKFITenergyvsdisplacement.py or CHECKFITforcevsdisplacement.py. If PLOTcheckderiv.py is executed after CHECKFITenergyvsdisplacement.py, it will give information on the energy derivatives, if it is executed after CHECKFITforcevsdisplacement.py, force derivatives are considered. The script PLOTcheckderiv.py is executed as follows.
$ PLOTcheckderiv.py YMIN YMAX
The previous plots can be used to determine the best range of displacements and order of polynomial fit. By analyzing the plot, we note that curves corresponding to the higher order of the polynomial used in the fit show a horizontal plateau at about 1584 cm^{1}. This can be assumed to be the converged value for the second derivative, from the point of view of the fit (see the analogue case in Energy vs. strain calculations).
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 6. for the explanation).
iii) PLOTstatus.py
Python visualization tool for the root mean squared (RMS) deviations of the effective SCF potential as a function of the iteration number during the SCF loop. It is executed as follows.
$ PLOTstatus.py LABEL
iv) PLOTatomforce.py
This script allows for the visualization of the forcevsdisplacement curve. It is executed as follows.
$ PLOTatomforce.py
6. Postprocessing: How to derive the optical phonon frequency at q=0 (Γ point)
i) Using the energyvsdisplacement curves
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 and a phonon pattern corresponding to the optical phonon at the Γ point, the displacements of the 2 atoms in the unit cell can be given as
u_{1}(R) = (0, 0, 0)
u_{2}(R) = a (u, u, u)
where a is the lattice constant. The same displacements hold in all unit cells of the crystal. Therefore, for an optical phonon at q=0, the total energy per unit cell of the crystal can be written as (see details here):
E({u}) = E_{0} + 3/4 m a^{2} (ω_{opt})^{2} u^{2}
Using appropriate values for m and a and their units, the phonon frequency of the optical mode at the Γ point for the diamond structure 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.
ii) Using the forcevsdisplacement curves
Using the same phonon displacement pattern as above, the component in the x direction of the force acting on the second atom is (within the harmonic approximation):
F_{2x}({u}) = 1/2 m a (ω_{opt})^{2} u
This force obviously vanishes in the equilibrium configuration (u=0). Similarly to the previous case, the frequency is obtained from the plateau value of the firstorder derivative at zero displacement of the forcevsdisplacement curves.
Exercises
 The frequency values which are obtained in this example using the quoted parameters, are not realistic. Check the convergence of the optical phonon frequency for diamond with respect of the value of the ngridk and rgkmax attributes using
 derivatives of the total energy,
 derivatives of the force acting on the atoms.
 Repeat all the calculation for Silicon.