Phonons At Γ In Diamond Structure Crystals¶
by Pasquale Pavone for exciting neon
(Jupyter notebook by Jennifer Peschel, Mara Voiculescu, & Martin Kuban)
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 energy-vs-displacement and force-vs-displacement curves at zero displacement.
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 tools/excitingjupyter/venv/excitingvenv/bin/activate
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 run_diamond_phonon_g.
%%bash
mkdir -p run_diamond_phonon_g
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
%%bash
cd run_diamond_phonon_g
python3 -m excitingscripts.setup.excitingroot
cd ..
Furthermore
do not perform any structural optimization.
perform the calculation of the force acting on the atoms by adding the
tforce="true"attribute to thegroundstateelement:
...
<groundstate
...
tforce="true">
</groundstate>
...
In order to generate input files for a series of different structures you have to run the script excitingscripts.setup.diamond_phonon_g. Notice that the script excitingscripts.setup.diamond_phonon_g 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 (in this example gamma_1).
If the directory name is not given explicitly, the script excitingscripts.setup.diamond_phonon_g generates a directory called by default workdir.
Very important: The working directory is overwritten each time you execute the script excitingscripts.setup.diamond_phonon_g. Therefore, choose different names for different calculations.
The following command line arguments are also required:
- The first argument (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 argument (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.
%%bash
cd run_diamond_phonon_g
python3 -m excitingscripts.setup.diamond_phonon_g 0.025 11 -w gamma_1
cd ..
2. Execute the Calculations¶
To execute the series of calculations with input files created by excitingscripts.setup.diamond_phonon_g you have to run the script excitingscripts.execute.diamond_phonon. If a name for the working directory has been specified, then you must give it here, too.
%%bash
cd run_diamond_phonon_g
python3 -m excitingscripts.execute.diamond_phonon -w gamma_1
cd ..
After the complete run, inside the directory gamma_1 the results of the calculations are contained in the subdirectories displ_*d* where d corresponds to the displacement value. The data for the energy-vs-displacement and force-vs-displacement curves are contained in the file phonon_results.json. More precisely, the force values correspond to the x component of the force acting on the second atom in the unit cell as a function of the displacement.
3. Post-Processing: Extract Energy Derivatives¶
In order to analyse the results of the calculations, you first have to move to the working directory. At this point, you can use the python script excitingscripts.checkfit.energy_vs_displacement for extracting the derivatives of the energy-vs-displacement curves at zero displacement.
The first argument (in our example 0.025) represents the absolute value of the maximum displacement for which we want to perform the analysis. The second argument (2) is the order of the derivative that we want to obtain. Finally, the third argument (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 the values of the frequency, i.e., combinations involving the square root of the derivative (see Section 6. for the explanation).
The script excitingscripts.checkfit.energy_vs_displacement generates the output file checkfit_energy_results.json, which can be used in the post-processing analysis.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.checkfit.energy_vs_displacement 0.025 2 12.01
cd ../..
The output should look should look something like this:
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.11 [cm-1]
Polynomial of order 3 ==> 1602.12 [cm-1]
Polynomial of order 4 ==> 1584.17 [cm-1]
Polynomial of order 5 ==> 1584.17 [cm-1]
Polynomial of order 6 ==> 1584.04 [cm-1]
Polynomial of order 7 ==> 1584.04 [cm-1]
Polynomial of order 2 ==> 48.0302 [THz]
Polynomial of order 3 ==> 48.0302 [THz]
Polynomial of order 4 ==> 47.4922 [THz]
Polynomial of order 5 ==> 47.4922 [THz]
Polynomial of order 6 ==> 47.4885 [THz]
Polynomial of order 7 ==> 47.4885 [THz]
4. Post-Processing: Extract Force Derivatives¶
The procedure for extracting the derivatives of the force-vs-displacement curves is very similar to the one exposed in the previous Section for the energy-vs-displacement curves. You can use the python script excitingscripts.checkfit.force_vs_displacement for extracting the derivatives of the force-vs-displacement curves at zero displacement. The meaning of the input entries for this script is the same as for excitingscripts.checkfit.force_vs_displacement (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.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.checkfit.force_vs_displacement 0.025 1 12.01
cd ../..
The output should look should look something like this:
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 [cm-1]
Polynomial of order 2 ==> 1605.96 [cm-1]
Polynomial of order 3 ==> 1580.60 [cm-1]
Polynomial of order 4 ==> 1580.60 [cm-1]
Polynomial of order 5 ==> 1580.62 [cm-1]
Polynomial of order 6 ==> 1580.62 [cm-1]
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. Post-Processing: Visualization Tools¶
All the scripts mentioned here must be executed in the directory where the phonon_results.json, checkfit_energy_results.json and checkfit_force_results.json files are located.
i) excitingscripts.plot.energy
This script allows for the visualization of the energy-vs-displacement curve. It is executed as follows.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.plot.energy
cd ../..
In the following, we display the result for the example mentioned above.
ii) excitingscripts.plot.checkderiv
This is a very important tool that allows to represent the dependence of the calculated derivatives of the energy-vs-displacement and force-vs-displacement 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").
If excitingscripts.plot.checkderiv is executed after excitingscripts.checkfit.energy_vs_displacement and by including the command line argument "energy", it will give information on the energy derivatives, if it is executed after excitingscripts.checkfit.force_vs_displacement and by including the command line argument "force", force derivatives are considered. The script excitingscripts.plot.checkderiv is executed as follows.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.plot.checkderiv energy
cd ../..
One or two optional entries can be specified on the calling line. Assigning numerical values to these two entries, you can set the minimum (YMIN) and the maximum (YMAX) value on the vertical axis, respectively. An example of the script output for the energy derivatives is the following.
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 the tutorial 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 energy-vs-displacement curves, but the values of the frequency, i.e., combinations involving the square root of the derivative (see Section 6. for the explanation).
iii) excitingscripts.plot.status
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.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.plot.status displ_-0.005
cd ../..
Here, the entry on the command line must be specified and represents the name of the directory where an exciting calculation has been or is still running. The excitingscripts.plot.status is particularly useful in the latter case because one can follow "live" the convergence of the calculation. An example of the output of the script is the following.
iii) excitingscripts.plot.atomforce
This script allows for the visualization of the force-vs-displacement curve. It is executed as follows.
%%bash
cd run_diamond_phonon_g/gamma_1
python3 -m excitingscripts.plot.atomforce
cd ../..
In the following, we display the result for the example mentioned above.
6. Post-Processing: How to Derive the Optical Phonon Frequency at q=0 (Γ point)¶
i) Using the energy-vs-displacement 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({\mathbf{u}}) = E_0 + \frac{1}{2} \sum_{\mathbf{R, R'}}\mathbf{u}(\mathbf{R})\: \mathbf{\Phi} (\mathbf{R, R'}) \:\mathbf{u}(\mathbf{R'}) + O(\mathbf{u}^3) $$where $E_0$ is the energy corresponding to the equilibrium structure, $\mathbf{u}(\mathbf{R})$ is an atomic displacement in the cell $\mathbf{R}$, $\mathbf{\Phi} (\mathbf{R, R'})$ is the inter-atomic force-constants matrix, and $O(\mathbf{u}^3)$ includes all orders higher than two in the displacements $\mathbf{u}$. For small displacements, terms beyond the harmonic approximation (which retains only terms that are quadratic in $\mathbf{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
$$ \mathbf{u}_1(\mathbf{R}) = (0, 0, 0)$$$$ \mathbf{u}_2(\mathbf{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({\mathbf{u}}) = E_0 + \frac{3}{4} m \: a^2 \: (\omega_{\text{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 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. The optical phonon frequency is obtained from the plateau value of the second-order derivative at zero displacement of the energy-vs-displacement curves. For more details on the procedure used for extracting numerical derivatives, see the tutorial Energy vs. strain calculations, where the same approach is used for calculating elastic constants.
i) Using the energy-vs-displacement 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}({\mathbf{u}}) = -\frac{1}{2} m \: a \: (\omega_{\text{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 first-order derivative at zero displacement of the force-vs-displacement 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
ngridkandrgkmaxattributes using- derivatives of the total energy,
- derivatives of the force acting on the atoms.
Repeat all the calculations for Silicon.