Lattice Dynamics of Diamond and Zincblende-Structure Crystals¶
by Sebastian Tillack & Pasquale Pavone for exciting sodium
(Jupyter notebook by Mara Voiculescu & Martin Kuban)
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 examples of diamond and cubic Boron nitride.
0. Before Starting¶
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:
source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate
In the following we use the conversion factor 1 Ha ≈ 2.194746 x 105 cm-1.
1. Theoretical Background¶
$\Rightarrow$ Theory of linear response for phonons
Within density-functional perturbation theory (DFPT) we calculate the linear response of the Kohn-Sham (KS) system upon an external perturbation $\delta$. We express the quantities of the perturbed system as an expansion in the perturbation operator and write down the KS equations for the perturbed system
\begin{equation} \tag{1} \left(\hat{\bf h}\phantom{}^0 + \lambda \,\delta \hat{\bf h} + \dots\right) \left(\psi_{n{\bf k}}^0({\bf r}) + \lambda\, \delta \psi_{n{\bf k}}({\bf r}) + \dots\right) = \left(\epsilon_{n{\bf k}}^0 + \lambda \,\delta \epsilon_{n{\bf k}} + \dots\right) \left(\psi_{n{\bf k}}^0({\bf r}) + \lambda \,\delta \psi_{n{\bf k}}({\bf r}) + \dots\right) \;. \end{equation}
Collecting terms linear in $\lambda$ we find the so called Sternheimer equation
\begin{equation} \tag{2} \left[\hat{\bf h}\phantom{}^0 - \epsilon_{n{\bf k}}^0\right] \,\delta \psi_{n{\bf k}}({\bf r}) = -\left[\delta \hat{\bf h} - \delta \epsilon_{n{\bf k}}\right] \,\psi_{n{\bf k}}^0({\bf r}) \;, \end{equation}
from which the linear response of the wavefunctions $\delta \psi_{n{\bf k}}$ can be obtained. Equivalently to the KS equations, the Sternheimer equation has to be solved self consistently, since the Hamiltonian response implicitly depends on the wave function response via the density and potential response:
\begin{equation} \tag{3} \delta \psi_{n{\bf k}} \rightarrow \delta n({\bf r}) \rightarrow \delta V({\bf r}) \rightarrow \delta \hat{\bf h} \rightarrow \delta \psi_{n{\bf k}} \;. \end{equation}
In the case of DFPT for phonons, the external perturbation is a collective phonon-like displacement of the nuclei from their equilibrium positions $\pmb{\tau}_{\kappa{\bf R}}$
\begin{equation} \tag{4} \delta = \delta^{\bf q}_{\kappa \alpha} \quad:\quad \pmb{\tau}_{\kappa {\bf R}} \longrightarrow \pmb{\tau}_{\kappa {\bf R}} + {\rm e}^{{\rm i}{\bf q}\cdot{\bf R}} \delta \pmb{\tau}_{\kappa} \;, \end{equation}
where $\kappa$ and ${\bf R}$ label atoms and unit-cells, respectively, ${\bf q}$ is the phonon wave vector and $\delta \pmb{\tau}_{\kappa}$ is an infinitesimal displacement vector.
In order to exploit symmetries, we do not displace every atom individually in the three Cartesian directions (canonical basis) but we use special symmetry adapted displacement patterns (irreducible representations, irreps):
\begin{equation} \tag{5} \delta^{\bf q}_{Id} = \sum_{\kappa,\alpha} p_{\kappa \alpha, Id}({\bf q})\, \delta^{\bf q}_{\kappa \alpha} \;, \end{equation}
where $I$ labels an irrep (for a given ${\bf q}$) and $d$ labels a so-called member of the irrep.
Once we have converged the Sternheimer equation for the density and potential response, we can use the results to compute in a single step the linear response of the atomic forces ${\bf F}_{\kappa}$ which directly correspond to the entries of the dynamical matrix
\begin{equation} \tag{6} \delta^{\bf q}_{\kappa \alpha} F_{\kappa' \beta} = \sqrt{M_\kappa\, M_{\kappa'}}\, D_{\kappa \alpha, \kappa' \beta}({\bf q}) \;. \end{equation}
2. Phonon Properties for Diamond¶
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.
%%bash
mkdir -p run_diamond_phonons
<input>
<title>Diamond: Phonons</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal scale="6.749908">
<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="4 4 4"
swidth="0.0001"
rgkmax="4.0"
gmaxvr="14"
tforce="true">
<libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
</groundstate>
<phonons
do="fromscratch"
method="dfpt"
ngridq="2 2 2"
polar="false">
<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>
The input block that describes a phonon calculation is specified by the element phonons. In exciting, there are two different methods available to compute phonons: the supercell method and DFPT. The latter is the default and what will cover in the first part of this tutorial. The dynamical matrices are calculated for a regular grid of phonon wavevectors q specified by the attribute ngridq inside phonons. 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. Afterwards, the phonon frequencies and eigenvectors at any wavevector q can be obtained by Fourier interpolation. Note that in the case of a DFPT calculation, the k- and q-grid need to be commensurate, i.e., ngridk must be an integer multiple of ngridq. In the case of a supercell calculation, both grids are completely uncorrelated.
A novelty in the groundstate element is that we replace the attribute xctype by the new element libxc. By that, we specify to take the exchange-correlation potential from the libxc library instead of the exciting internal implementation. The reason for this change is that a DFPT calculation requires the xc-kernel (the second derivatives of the xc-energy) which is not natively implemented in exciting and is always obtained from the libxc library. For full consistency, we recommend to also use the xc-potential (first derivative of the xc-energy) from libxc. However, in principle you can also use the exciting internal implementations of the xc-potential using the attribute xctype. In this case, only the xc-kernel of the corresponding potential is obtained from libxc.
Further attributes of the element phonons are the following.
method, that controls the method employed to compute phonons (default is"dfpt").reduceq, which uses crystal symmetry to reduce the q-points (default is"true").polar, that specifies whether a polar or nonpolar material is to be calculated (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.3.
After you have created the input file input.xml, you can run exciting by using the commands
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.setup.excitingroot
python3 -m excitingscripts.execute.single
cd ..
The calculation is completed if the file PHONON.OUT has been created.
2.2 Output Files of a Phonon Calculation
The main result of a phonon calculation are the dynamical matrices. They are given in the files DYN_fileext.OUT, where fileext means an extension of the filename having the structure Q############_S##_A###_P#.OUT and specifying the q vector on the grid (Q…) species and index of the displaced atom (S… and A…) and polarization direction (P…). Additionally, a number of subdirectories is created. In a DFPT calculation, their names are Q############_I###, where the first part specifies the q-point in the grid and the second part labels the irrep at that specific q-point. Irreps are special symmetry adapted displacement patterns. Inside of each subdirectory, you will find the following files:
| filename | description |
|---|---|
| PHONON_INFO_fileext.OUT | General information about the run. Similar to INFO.OUT for the groundstate. |
| DYN_fileext_#.OUT | Force response (or dynamical matrix row) for the #-th member of the irrep. |
| PHONON_DRHO_fileext_#.OUT | Electron density response for the #-th member of the irrep. |
| PHONON_DVEFF_fileext_#.OUT | Effective potential response for the #-th member of the irrep. |
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.
2.2.1 PHONON_INFO_fileext.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.
k-points are reduced by symmetries.
Response quantities are symmetrized.
q-points and symmetries:
iq q-vector (lattice coordinates) N_sym N_k N_irrep (size)
--------------------------------------------------------------------------------
1 ( 0.000000, 0.000000, 0.000000) 48 8 2 (3*,3)
2 ( 0.500000, 0.000000, 0.000000) 12 13 4 (2,1,1,2)
3 ( 0.500000, 0.500000, 0.000000) 16 13 3 (2,2,2)
current task / total (remaining) number of tasks: 1/ 9 (9)
load of current task / total (remaining) load: 24/ 204 (204)
procs for current task / total number of procs: 1/ 1
average / maximum load per proc: 204/ 204
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) q-points 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 k-points in the irreducible Brillouin zone obtained with the symmetries in the small group. N_irrep denotes the number of irreps for the given q-point and the values in parenthesis give their respective size (or equivalently the number of members per irrep). The q-point 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 5.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 self-consistency cycle of the Sternheimer equation and shows the change in the potential response for each iteration:
********************************************************************************
* scf loop for Sternheimer equation started *
********************************************************************************
Sternheimer SCF loop iteration 1
maximum change in potential response: 0.960508E-01
Fermi energy response: 0.00000
Time spent for iteration (seconds): 1.20
Sternheimer SCF loop iteration 2
maximum change in potential response: 0.646932E-01
Fermi energy response: 0.00000
Time spent for iteration (seconds): 3.40
...
2.2.2 PHONON.OUT
If you have specified a qpointset in the input file, 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 excitingscripts.convert.au2invcm to obtain frequencies expressed in units of cm-1.
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.convert.au2invcm
cd ..
- 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.3 Properties from the Phonon Calculation
To obtain other phonon properties, additional subelements can be inserted inside the element phonons, as explained in the following.
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 ngridqint="70 70 70" nwdos="1000" nsmdos="2"/>
...
Here, phonon frequencies are computed on a dense grid (ngridqint) and collected to obtain the DOS, which is written to the file PHDOS.OUT. The number of steps between the lowest and highest frequency (N$_{\omega}$) is given by nwdos, nsmdos is the number of times the data are smoothened.
The modified input file in this case is shown below.
<input>
<title>Diamond: Phonons</title>
<structure speciespath="$EXCITINGROOT/species">
...
</structure>
<groundstate
do="skip"
.../>
<phonons
do="skip"
...>
<phonondos
ngridqint="70 70 70"
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 2.1) and to work in the same directory. If this is the case, in order to perform the calculation, you type
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.execute.single
cd ..
The data in PHDOS.OUT can be plotted directly. For this, you can use the script excitingscripts.plot.dos, which is fully described in The python script "plot.dos". Execute the script as follows.
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.plot.dos -p -e 0 1400
cd ..
The script produces a PNG (PLOT.png) output, which should be similar to the one displayed here.
A better converged result (using an 8×8×8 k-grid, rgkmax = "7.0", and 4 q-points in each direction) can be found below.
$\Rightarrow$ Show converged plot
Hint: Run the converged calculation in a different directory, e.g., /home/tutorial/diamond-phonons/conv. Copy here the previous input.xml and modify the computational parameters. The result should look like the following.
Thermodynamic properties
The zero-point energy, $E_{zp}$, and further thermodynamic properties can be computed from the phonon DOS $g(\omega)$ (see Input Reference for their definition for further details).
- the vibrational free energy $F_{vib}= U_{vib} - TS_{vib}$,
- the vibrational internal energy $U_{vib}$,
- the entropic contribution to the vibrational free energy $TS_{vib}$,
- the vibrational entropy $S_{vib}$
- the heat capacity $C_v$.
In order to calculate and visualize these thermodynamic properties, you can use the script excitingscripts.free_from_dos as follows
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.free_from_dos 0 1400 160 eV
cd ..
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 generates the output files F_vib, U_vib, TS_vib, S_vib, C_v, which can be visualized with standard tools. Examples using the script excitingscripts.plot.files are given in the following:
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.plot.files -f F_vib U_vib TS_vib -lx 'Temperature [K]' -ly 'Free energy [eV]' -g -mtx 4 -s 0.65 1.3 -pn 'PLOT-free-energies'
python3 -m excitingscripts.plot.files -f S_vib -lx 'Temperature [K]' -ly 'Entropy [meV/K]' -g -pn 'PLOT-entropy'
python3 -m excitingscripts.plot.files -f C_v -lx 'Temperature [K]' -ly 'Heat capacity [$k_B$]' -g -pn 'PLOT-heat-capacity'
cd ..
$\Rightarrow$ Show converged plots
Zero-point energy is 3.5905e-01 [eV]
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.
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.execute.single
cd ..
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 excitingscripts.plot.band_structure (see The python script "plot.band_structure"), can be used to produce a plot of the phonon-dispersion relations.
%%bash
cd run_diamond_phonons
python3 -m excitingscripts.plot.band_structure -p -e 0 1400
cd ..
$\Rightarrow$ Show converged plot
Using the path Γ-(K)-X-Γ-L this result is obtained.
$\Rightarrow$ Show converged plot with standard path
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). 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>
...
The standard path is visualized in the following figure:
Using these settings, the following plot is obtained.