by Pasquale Pavone & Stefan Kontur for exciting oxygen
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 phonon-dispersion relations. You will work through the example of diamond.
|
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.
- FREE-fromdos.py: Python script for calculating thermodynamic properties.
- PLOT-dos.py: Python visualization tool for the phonon density of states.
- PLOT-band-structure.py: Python visualization tool for phonon-dispersion curves.
- PLOT-files.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
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/diamond-phonons.
$ cd /home/tutorials
$ mkdir diamond-phonons
$ cd diamond-phonons
1. Set up an input file for a phonon calculation
If you already completed Phonons at Γ in diamond-structure crystals you can copy the input file into diamond-phonons as a starting point and modify it as discussed below.
1.1) Define structural parameters
First, the parameter defining your structure should be defined. In this example, we consider carbon in the diamond phase. The structure part of the input file looks then like the following.
... <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> ...
The chosen lattice parameter of 6.7468 refers to the converged equilibirum parameter for diamond, as obtained with the PBE exchange-correlation functional.
1.2) Parameters for the ground-state calculation
Then, parameters for performing the SCF have to be specified. In this example, we use the GGA functional of Perdew, Burke, and Ernzerhof (PBE), other choices can be found in the description of the attribute xctype of the element groundstate.
... <groundstate do="fromscratch" ngridk="2 2 2" swidth="0.0001" rgkmax="4.0" gmaxvr="14" maxscl="25" xctype="GGA_PBE"/> ...
Note that this example of input contains values for the parameters ngridk and rgkmax which allow for fast calculations but should not be considered realistic. A detailed description of how to achieve convergence with respect to these parameters can be found in Simple convergence tests.
1.3) Set up the supercell phonon calculation
Simple frozen-phonon calculations can be performed only at high-symmetry points in the Brillouin zone (see, e.g., Phonons at Γ in diamond-structure crystals and Phonons at X in diamond-structure crystals). With exciting, we can also perform phonon calculations at any point in the Brillouin zone using the supercell method and phonon interpolation. The input block in which this method is used is the one specified by the element phonons. In supercell phonon calculations, the size of the supercell is specifyed by the periodicity (the wavevector q) of the phonon displacements which are considered. For instance, as mentioned in Phonons at X in diamond-structure crystals, the supercell corresponding to the X point has a size which is twice the one of the unperturbed crystal.
In general, for calculating phonon properties, you must specify the attribute ngridq inside phonons. Setting ngridq = "n1 n2 n3", we are asking exciting to calculate all the phonons corresponding to all the q-points belonging to a "n1 n2 n3" mesh in reciprocal space. In this case, exciting performs calculation for supercell defined by multiplying the crystal basis vectors of the unit cell by n1, n2, and n3, respectively. The q-points mesh is defined in the same way as the k-points grid (specified by ngridk), which is used in calculating the electronic ground state. However, the choices of the k-points and q-points meshes are completely uncorrelated.
In the following example, we choose ngridq = "2 2 2", which corresponds to a supercell with all crystal basis vectors doubled with respect to the unit cell of the unperturbed crystal.
... <phonons do="fromscratch" ngridq="2 2 2"> ... </phonons> ...
Further attributes of the element phonons are the following.
- deltaph, that controls the step size for atom displacements (default is 0.03 Bohr).
- reduceq, which uses crystal symmetry to reduce the q-points (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 q-points specified in each subelement qpoint. In the example above, calculations will be performed for Γ, X, and L phonons (in the standard notation for high-symmetry points in the Brillouin zone of fcc crystals). Note that q-points are given in reciprocal lattice coordinates. Additional properties are discussed in Section 2.
1.4) The basic input file for a phonon calculation
We are now ready to run the actual phonon calculation. Create a subdirectory c-test and move inside it to run the calculations.
$ mkdir c-test
$ cd c-test
Jobs which runs a phonon calculation are in general a bit lengthy, because they use supercells which are larger than the primitive cell of the crystal. The present calculation should take about 15 minutes on a recent CPU with the loose parameters specified. Note that they can also be run in parallel. The complete input file for this calculation is given in the following.
<input> <title>Diamond: 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" xctype="GGA_PBE"/> <phonons do="fromscratch" 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>
Please, remember that after creating the input file input.xml, the path for the species file must be set using the command
$ SETUP-excitingroot.sh
Now, you can run exciting by using the command
$ 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
1.5) Output files of a phonon calculation
The dynamical matrix is set up from a series of supercell ground-state calculations. In a typical supercell calculation of this series, all atoms which are equivalent to a given atom in the primitive cell of the unpertubed crystal are displaced along a cartesian direction. 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).
After a successful phonon calculation, a number of subdirectories are created. The names of these directory 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.
A typical phonon calculation produces in addition the following files.
|
Here, fileext means an extension of the filename, again naming the corresponding q-point of the grid, species, atom, and polarization.
As indicated above, the frequencies of the phonon modes at the required q-points 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 CONVERT-au2invcm.py to obtain frequencies expressed in units of cm-1.
$ CONVERT-au2invcm.py
############################################
--------------------------------------------
Total number of atoms : 2
Total number of q-points : 3
--------------------------------------------
q-point 1 is 0.0 0.0 0.0
mode 1 frequency: 0.0000 cm-1
mode 2 frequency: 0.0000 cm-1
mode 3 frequency: 0.0000 cm-1
mode 4 frequency: 1580.5232 cm-1
mode 5 frequency: 1580.5232 cm-1
mode 6 frequency: 1580.5232 cm-1
q-point 2 is 0.5 0.5 0.0
mode 1 frequency: 523.8799 cm-1
mode 2 frequency: 523.8799 cm-1
mode 3 frequency: 1096.3671 cm-1
mode 4 frequency: 1096.3671 cm-1
mode 5 frequency: 1128.0149 cm-1
mode 6 frequency: 1128.0149 cm-1
q-point 3 is 0.5 0.5 0.5
mode 1 frequency: 459.4405 cm-1
mode 2 frequency: 459.4405 cm-1
mode 3 frequency: 869.2564 cm-1
mode 4 frequency: 1238.7319 cm-1
mode 5 frequency: 1400.8122 cm-1
mode 6 frequency: 1400.8122 cm-1
############################################
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. 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.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 ngrdos="300" nwdos="1000" nsmdos="2"/> ...
Here, phonon frequencies are computed on a dense grid (ngrdos
The complete input file in this case is shown below.
<input> <title>Diamond: PhononDOS</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="skip" ngridk="2 2 2" swidth="0.0001" rgkmax="4.0" gmaxvr="14" xctype="GGA_PBE"/> <phonons do="skip" ngridq="2 2 2"> <phonondos ngrdos="300" 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 1.4) and to work in the same directory. If this is the case, in order to perform the calculation, you type
$ SETUP-excitingroot.sh
$ time exciting_smp &
Data in PHDOS.OUT can be plotted directly. You can use the script PLOT-dos.py (see The python script "PLOT-dos.py"), by typing on the command line
$ PLOT-dos.py -p
The script produces a PNG (PLOT.png) output, which should be similar to the one displayed here.
A better converged result (using an 8
Hint: Run converged calculation in a different directory, e.g., /home/tutorial/diamond-phonons/c-conv.
$ mkdir /home/tutorial/diamond-phonons/c-conv
$ cd /home/tutorial/diamond-phonons/c-conv
Copy here the previous input.xml and modify the computational parameters. The result should look like the following.
2.2) Thermodynamic properties
The zero-point energy, Ezp, and further thermodynamic properties can be computed from the phonon DOS g
- the vibrational free energy Fvib= Uvib -TSvib,
- the vibrational internal energy Uvib,
- the entropic contribution to the vibrational free energy TSvib,
- the vibrational entropy Svib
- the heat capacity C,,v.
In order to calculate and visualize these thermodynamic properties, you can use the script FREE-fromdos.py as follows
$ FREE-fromdos.py 0 1600 160 eV
Zero-point energy is 3.5041e-01 [eV]
Created PNG output "PLOT-free-energies.png"
Created PNG output "PLOT-entropy.png"
Created PNG output "PLOT-heat-capacity.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 zero-point energy is written on the screen output. The script also generates three PNG output files: PLOT-free-energies.png, PLOT-entropy.png, and PLOT-heat-capacity.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.
Zero-point energy is 3.5850e-01 [eV]
2.3) Phonon-dispersion relations
Phonon-dispersion 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 PLOT-band-structure.py (see The python script "PLOT-band-structure.py"), can be used to produce a plot of the phonon-dispersion relations.
$ PLOT-band-structure.py -p -a KS
The resulting PNG file PLOT.png looks like the following.
Converged results as mentioned above can be seen below.
Using the path Γ-(K)-X-Γ-L this result is obtained.
A complete visualization of the phonon dispersion of diamond can be achieved using the standard path in reciprocal space for fcc crystals introduced by Setyawan and Curtarolo, Comp. Mat. Sci. 49, 299-312 (2010), which is visualized in the following figure:
The coordinates of the high-symmetry points of the first Brillouin zone belonging to this path are given below.
... <phonondispplot> <plot1d> <path steps="400"> <point coord="0.000000 0.000000 0.000000" label="Gamma"/> <point coord="0.500000 0.000000 0.500000" label="X"/> <point coord="0.500000 0.250000 0.750000" label="W"/> <point coord="0.375000 0.375000 0.750000" label="K"/> <point coord="0.000000 0.000000 0.000000" label="Gamma"/> <point coord="0.500000 0.500000 0.500000" label="L"/> <point coord="0.625000 0.250000 0.625000" label="U"/> <point coord="0.500000 0.250000 0.750000" label="W"/> <point coord="0.500000 0.500000 0.500000" label="L"/> <point coord="0.375000 0.375000 0.750000" label="K"/> </path> </plot1d> </phonondispplot> ...
Using these settings, the following plot is obtained.
2.4) Other useful features
- The element interpolate triggers an interpolation on a given q-point grid. For every point in the grid, which may be offset from the origin, phonon frequencies and (optionally) eigenvectors can be obtained, and written to the file PHONON_INTERPOLATE.OUT.
- The element reformatdynmat (without any attributes) requests an output of the dynamical matrix in the form of a 3
× 3 block for each combination of atoms and each q-point (file DYNMAT.OUT). - The file DYNMAT_SUMRULE.OUT contains the dynamical matrix (in 3
× 3 blocks as above) after applying a correction, such that the matrix for q=0 satisfies the acoustic sum rule. This ensures that the q=0 dynamical matrix has 3 zero eigenvalues. The files DYNMAT_SYM.OUT and DYNMAT_SYM_SUMRULE.OUT contain the symmetrized matrices.
Exercise: Silicon in diamond structure
- The calculation presented here for diamond can be repeated for silicon. Use the PBE exchange-correlation functional, an equilibrium lattice parameter of 10.34 Bohr, and low values for both k- and q-point grids to accelerate the runs, as above.
- Compare your results with the "converged" one displayed here.
http://exciting.wikidot.com/oxygen-phonon-properties-diamond
Hello, my name is Dr Abderrahmane Reggad. I'm a 53 year old and I am a researcher in materials' sciences. 
0 Comments