for helium version of Exciting
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 energy-vs-displacement curves at zero displacement.
|
Table of Contents
|
1. Define relevant shell variables and download scripts
Very important: Before starting, the following shell variables must be set by the user:
EXCITINGROOT = Directory where exciting has been downloaded, e.g.: /home/user/exciting .
EXCITINGRUNDIR = Directory where the user runs all the calculations (must be existing!!), e.g.: /home/user/tempdir .
EXCITINGSCRIPTS = Directory where the scripts are downloaded (not necessary if the downloaded files are copied in /usr/bin), e.g.: /home/user/scriptsdir .
The setting of these variables can be done in a bash shell by typing (from now on the symbol $ will indicate the shell prompt):
$ export EXCITINGROOT=/local_path_to_exciting_root
$ export EXCITINGRUNDIR=/local_run_directory
$ export EXCITINGSCRIPTS=/local_tutorial_scripts_bin_directory
$ export PATH=$PATH:$EXCITINGSCRIPTS
The directory names /local_path_to_exciting_root, /local_run_directory, and /local_tutorial_scripts_bin_directory are dummy names which must be explicitly changed by the user to the appropriate ones.
As a second step, it is necessary to download the scripts which are used in this tutorial. Here is a list of these scripts with a short description, click on the script name to download it:
- SETUP-phonon-X.py: Python script for generating structures with displaced atoms.
- EXECUTE-phonon-X.sh: (Bash) shell script for running a series of exciting calculations.
- CHECKFIT-energy-vs-displacement.py: Python script for extracting the derivatives at zero displacement of the energy-vs-displacement curves.
- GNU-energy: Gnuplot visualization tool for the energy-vs-displacement curves.
- GNU-status: Gnuplot visualization tool for the RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
- GNU-check: Gnuplot visualization tool for the fit of the derivatives of the energy-vs-displacement curves at zero displacement.
All the downloaded scrips must be executable, if this is not the case use, e.g.:
$ chmod u+x SETUP-phonon-X.py
Requirements: Bash shell. Python numpy, lxml, and sys libraries. Gnuplot (visualization tools).
2. 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 energy-vs-displacement curves for carbon in the diamond structure. Thus, we will create a directory diamond and we move inside it:
$ mkdir diamond
$ cd diamond
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:
<?xml version="1.0" encoding="UTF-8"?> <?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?> <input xsi:noNamespaceSchemaLocation="/local_path_to_exciting_root/xml/excitinginput.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsltpath="/local_path_to_exciting_root/xml/"> <title>Diamond: equilibrium structure</title> <structure speciespath="/local_path_to_exciting_root/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="4 4 3" xctype="GGAPerdew-Burke-Ernzerhof" stype= "Methfessel-Paxton 1" swidth="0.001" gmaxvr="18" lmaxvr="8" rgkmax="7.0" lradstep="2"> </groundstate> </input>
Please, remind that the input file for an exciting calculation must be always called input.xml.
Be sure to:
- set the correct path for the species directory (/local_path_to_exciting_root) to the one pointing to the place where the exciting directory is placed:
... <input xsi:noNamespaceSchemaLocation="/local_path_to_exciting_root/xml/excitinginput.xsd" ... xsltpath="/local_path_to_exciting_root/xml/"> ... <structure speciespath="/local_path_to_exciting_root/species"> ...
- do not perform any structural optimization.
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, 1, 0) a / √2;
- (1, 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 crystal and Cartesian coordinates coincide):
- (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 SETUP-phonon-X 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 looks 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 should be also changed according to the new geometry. In particular, due to the fact that the crystal axis have now lengths which relate to each other like "1 1 √2", the values in ngridk should be choosen to relate each other as close as possible (notice that the values must be integers) to "1 1 1/√2". See the tutorial Eletronic-Structure Calculations for more details on this topic. Finally, the k-point mesh should be chosen in the original input file (and imported to the new ones) as, e.g.:
... <groundstate ngridk="4 4 3" ...
Other possible choices for the k-point grid are the following:
- "6 6 4"
- "8 8 6"
- "10 10 7"
- …
iii) Generation of input files for the different structures
In order to generate input files for a series of different structure you have to run the script SETUP-phonon-X, which will produce the following output on the screen:
$ SETUP-phonon-X.py
Enter maximum displacement umax [u/(lat. parameter)] >>>> 0.03
Enter the number of displacements in [-umax,umax] >>>> 13
------------------------------------------------------------------------
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.03) 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 (13) 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.
After running the script, a directory called workdir is created, which contains the input files corresponding to the different displacements.
Notice that the directory workdir will be overwritten each time you execute the script SETUP-phonon-X.py. Therefore, before setting up a new run please rename the workdir directory to a different name, e.g.:
$ mv workdir results-label
3. Execute the calculations
To execute the series of calculation with input files created by SETUP-phonon-X.py you have to run the script EXECUTE-phonon-X.sh:
$ EXECUTE-phonon-X.sh
Running exciting for file input-07.xml ----------------------------------
...
Run completed for file input-13.xml -------------------------------------
$
Notice that, due to the symmetry of the phonon pattern at X, the script executes only calculations with input files corresponding to positive (or vanishing) displacements. This is the reason why the first executed calculation is for the file input-07.xml instead of for input-01.xml.
After the complete run, the results of the calculation for the input file workdir/input.xml-i are contained in the subdirectory workdir/rundir-i where i is running from 01 (07 in this case) to the total number of strain values. The data for the energy-vs-displacement curves are contained in the files workdir/energy-vs-displacement.
4. Post-processing: Extract energy derivatives
In order to analyze the results of the calculations, you first have to move to the working directory.
$ cd workdir
At this point, you can use the python script CHECKFIT-energy-vs-displacement.py for extracting the derivatives of the energy-vs-displacement curves at zero displacement:
$ CHECKFIT-energy-vs-displacement.py
Enter maximum displacement for the interpolation >>>> 0.03
Enter the order of derivative >>>> 2
Enter atomic mass in [amu] >>>> 12.01
_____________________________________________
X-phonon-calculation
#############################################
Fit data------------------------------------
Maximum value of the displacement ==> 0.030
Number of displacement values used ==> 13
Fit results for the derivative of order 2
Polynomial of order 2 ==> 778.80 [cm-1]
Polynomial of order 3 ==> 778.80 [cm-1]
Polynomial of order 4 ==> 774.90 [cm-1]
Polynomial of order 5 ==> 774.90 [cm-1]
Polynomial of order 6 ==> 774.12 [cm-1]
Polynomial of order 7 ==> 774.12 [cm-1]
Polynomial of order 2 ==> 23.3479 [THz]
Polynomial of order 3 ==> 23.3479 [THz]
Polynomial of order 4 ==> 23.2309 [THz]
Polynomial of order 5 ==> 23.2309 [THz]
Polynomial of order 6 ==> 23.2076 [THz]
Polynomial of order 7 ==> 23.2076 [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.03) represents the absolute value of the maximum displacement for which we want to perform the calculation. 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 energy-vs-displacement curves, but a combination involving the square root of such quantities (see Section 6. for the explanation).
The script generates the output files check-energy-derivatives and order-of-derivative, which can be used in the post-processing analysis.
5. Post-processing: Visualization tools
All the scripts mentioned here must be executed in the directory where the energy-vs-displacement, check-energy-derivatives, and order-of-derivative files are located. The PostScript output file is always named gnu.ps.
i) GNU-energy
This script allows for the visualization of the energy-vs-displacement curve. It is executed as follows:
$ GNU-energy
The script generates the PostScript file gnu.ps. In the following, we display the result for the example mentioned above (carbon in the diamond structure, displacements up to 3% in units of the cubic lattice constant):
ii) GNU-check
This is a very important script that allows to represent the dependence of the calculated derivatives of the energy-vs-displacement curve on
- the range of points included in the fitting procedure,
- the maximum degree of the polynomial used in the fitting procedure.
The script is executed as follows:
$ GNU-check
An example of the script output is the following:
The most appropriate value for the derivative in this example is represented by the horizontal plateau (at about 774 cm-1) of the curves corresponding to the highest order fit.
NOTICE that, in this example, the values which are plotted in the picture above are not directly the actual second derivative of the energy-vs-displacement curves, but a combination involving the square root of such quantities (see Section 6. for the explanation).
iii) GNU-status
Gnuplot visualization tool for the RMS deviations of the effective SCF potential as a function of the iteration number during the SCF loop. It is executed as follows:
$ GNU-status
Enter label [r or 01,02,...] >>>> 13
$
In this example, the input entry, preceded by the symbol »» (in our example 01) represents the label of the calculation for which you would like to plot the RMS deviations. In particular, the choice r refers to the currently running calculation and, e.g., 13 to the calculation already saved in the directory workdir/rundir-13. An example of the script output is the following:
Remark on the visualization scripts
All the visualization scripts accept two on-the-command-line entries specifying the y-axis limiting values. So that, e.g., the command
$ GNU-check 700 800
will visualize the results for the derivative in the y-axis range between 700 and 800 cm-1 .
6. Post-processing: 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
E({u}) = E0 + 1/2 ΣR,R' u(R) Φ(R,R') u(R') + O(u3)
where E0 is the energy corresponding to the equilibrium structure, u(R) is an atomic displacement in the cell R, and Φ(R,R') is the inter-atomic force-constants matrix. 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
- u1 = a (0, u, 0)
- u2 = a (0, u, 0)
- u3 = a (0, -u, 0)
- u4 = a (0, -u, 0)
LA phonon
- u1 = a (0, 0, u)
- u2 = a (0, 0, u)
- u3 = a (0, 0, -u)
- u4 = a (0, 0, -u)
TO phonon
- u1 = a (0, u, 0)
- u2 = a (0, -u, 0)
- u3 = a (0, -u, 0)
- u4 = a (0, u, 0)
LO phonon
- u1 = a (0, 0, u)
- u2 = a (0, 0, -u)
- u3 = a (0, 0, -u)
- u4 = 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:
E({u}) = E0 + m a2 ω(X)2 u2
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 energy-vs-displacement curves. This is the output of the fit procedure exposed above if the order of the derivative is taken equal to 2. One obtains, e.g., in our case for a TA mode
ωTA(X) = 774.12 cm-1 = 23.2076 THz
if the order of the fitting polynomial used for the fit is 4 and 13 different displacements equally spaced between -0.03 and +0.03 are considered. If the (more realistic) plateau value is used, as explained above, the value
ωTA(X) = 774 cm-1 = 23.20 THz
is obtained.