Exercise 2: Running Phonopy calculations via FHIvibes
Estimated total CPU time: 1 min
Warning
In the following exercises, computational settings including the reciprocal space grid (tag k_grid
), the basis set, and supercell's size, have been chosen to allow for a rapid computation of the exercises in the limited time and within the CPU resources available during the tutorial session. Without loss of generality, these settings allow to demonstrate trends of the lattice dynamics of materials.
In the production calculation, all computational parameters should be converged.
In this exercise, you will learn:
 How to use FHIvibes to set up and perform phonons calculation.
 How to compute phonon band structures using Phonopy.
 How to compute resulting properties like density of states (DOS), vibrational free energies and heat capacities.
 How to visualize phonon modes (optional, but highly recommended).
As input structure use the final relaxed structure from the previous exercise.
We start by generating the input file with general calculator settings:
vibes template calculator aims >> phonopy.in
In order to add settings for Phonopy just type:
vibes template phonopy >> phonopy.in
Section phonopy
We finally come to the section where we define the parameters for the Phonopy calculation.
cat phonopy.in
...
[phonopy]
supercell_matrix: [1, 1, 1]
displacement: 0.01
is_diagonal: False
is_plusminus: auto
symprec: 1e05
q_mesh: [45, 45, 45]
workdir: phonopy
The ransformation matrix between primitive cell
and supercell is provided in the supercell_matrix
variable.
One could also write it as [1, 0, 0, 0, 1, 0, 0, 0, 1] or
[[1, 0, 0], [0, 1, 0], [0, 0, 1]], which would be equivalent for FHIvibes parser, which
understands any array of numbers that can be converted to
3\(\times\)3 shape.
The supercell matrix \(M_\text{s}\) will be used to generate the lattice of the supercell from the lattice of the primitive unit cell by matrix multiplication ^{1}:
Here, \(\boldsymbol{a}_\text{u}\), \(\boldsymbol{b}_\text{u}\), \(\boldsymbol{c}_\text{u}\) are the lattice vectors of the (primitive) unit cell and \(\boldsymbol{a}_\text{s}\), \(\boldsymbol{b}_\text{s}\), \(\boldsymbol{c}_\text{s}\) label the lattice vectors of the supercell respectively.
For the start, we do not change the trivial supercell matrix that corresponds to the unit matrix. In this case, the "supercell" is just the primitive unit cell. We will investigate the influence of the supercell size and shape later on.
In Phonopy ^{8}, force constants are generated based on the finite displacement method ^{6}. Crystal symmetry is used to reduce the calculation cost and numerical noise of the force constants. Firstly, a symmetry reduced set of atomic displacements is generated. Then forces are computed for each structure from that set. After the atomic force calculations, the set of atomic displacements are expanded using the symmetry and then all the elements of force constans between atoms in a primitive cell and the supercell are fit to the symmetry expanded forces of atoms in supercells using Moore–Penrose pseudoinverse. This procedure may considered as a variant of ParlinskiLiKawazoe method ^{7}.
Run phonopy calculation
You may now perform calculations inside this folder by running
vibes run phonopy
The calculation should take only a few seconds and you should see a
folder phonopy
that contains your calculations, already
collected into the file called trajectory.son
. This file
contains all the necessary information from the DFT calculation, i.e.,
the structures with displaced atoms, and the according DFT forces. With
this information, FHIvibes can "feed" Phonopy and determine the
force constants and the dynamical matrices at arbitrary q points, from
which the phonon frequencies can be computed.
To postprocess the calculation, type:
vibes output phonopy phonopy/trajectory.son v
You should then obtain information about the phonon frequencies:
...
Frequencies at Gamma point:
q = [0. 0. 0.] (weight= 1)
# Mode Frequency
1 0.0000004 THz
2 0.0000004 THz
3 0.0000002 THz
4 15.4878848 THz
5 15.4878852 THz
6 15.4878857 THz
Congratulations, you have just performed your first phonon calculation and computed phonon frequencies at the \(\Gamma\) point (\(\mathbf{q} = 0\))^{2}! But why are there six frequencies and why do the first three are (numerically) equal to zero? Note that by default once the postprocessing of the output is done FHIvibes prints only the information about frequencies at the \(\Gamma\) point, whereas the way to obtain full spectrum would be discussed below.
Postprocessing, Band Structure, DOS etc.
We will now perform slightly more postprocessing of the Phonopy calculations. Before we progress futher, a mini recap on the periodic solids:
Brillouin Zone
We have just learned how to compute the Hessian of a given periodic system, and how to use this information to compute the dynamical matrix and its eigenvalues for the \(\Gamma\) point. As you know, the periodicity of the atomic positions in a lattice is reflected by the fact that it is sufficient to look at \(\mathbf{q}\) values within the first Brillouin zone of the reciprocal lattice: A wave vector \(\mathbf{q}\) that lies outside the first Brillouin zone corresponds to a wave whose wavelength is shorter than the distance between periodic images of the same atom. It can thus be represented equally well by a wave with longer wavelength, i. e. a smaller wave vector \(\tilde{{\mathbf{q}}}\) taken from the first Brillouin zone. ^{3}
High Symmetry Points
The Brillouin zone of our Silicon fcc diamond structure is displayed in the above Figure. The labelled points correspond to \(\mathbf{q}\) values of high symmetry. This means that there are symmetry operations in the point group of the lattice that leave this point invariant (up to a reciprocal space vector) ^{5}.
You can list the high symmetry points of the lattice of your geometry with FHIvibes by typing
vibes info geometry geometry.in vv
Please verify that you obtain the high symmetry points
...
Special k points:
G: [0. 0. 0.]
K: [0.375 0.375 0.75 ]
L: [0.5 0.5 0.5]
U: [0.625 0.25 0.625]
W: [0.5 0.25 0.75]
X: [0.5 0. 0.5]
Note that the list of points is given in fractional coordinates as coefficients of the reciprocal lattice. For the meaning of the Symbols \(\Gamma\), \(X\), etc., you can take a look at en.wikipedia.org/wiki/Brillouin_zone.
Band Structure
If we connect two or more \(\mathbf{q}\) points from the Brillouin zone, solve the eigenvalue problem for any \(\mathbf{q}\) point in between, and plot the obtained dispersions \(\omega (q)\) versus \(q\), we obtain the socalled phonon band structure. The band structure is typically computed for a path in the Brillouin zone that connects several or all of the high symmetry points. We will use the \(\mathbf{q}\) paths proposed by Setyawan and Curtarolo ^{3}, as shown above.
In order to compute the phonon band structure, we perform the following steps:

Run the Phonopy calculations with FHIvibes as you have just learned it (should be finished from the first part of this exercise).

Use
vibes output phonopy phonopy/trajectory.son bs
to perform the postprocessing including the band structure.
The last command will inform you that the phonons' spectrum was computed,
written, and plotted to files you fill find in phonopy/output
:
...
Extract phonopy results:
.. q_mesh: [45, 45, 45]
.. write primitive cell
.. write supercell
.. write force constants
.. plot band structure
.. write band structure yaml file
.. all files written to phonopy/output in 1.233s
Inspect the phonons spectrum in a pdf viewer before you proceed. Despite using a minimal "supercell", important features should already be visible. It should like the following figure.
DOS and Thermal Properties
After you managed to compute the band structure, we proceed to evaluate Eq. (5), i.e., how to compute the density of states and plot it. As before, you just need to perform the respective postprocessing, this time with
vibes output phonopy phonopy/trajectory.son dos
Phonopy evaluates the frequencies on a grid of
45\(\times\)45\(\times\)45 \(\mathbf{q}\)
points per default and uses the socalled Tetrahedron method to
interpolate between the points ^{4}. Afterwards it literally
counts the number of frequencies in bins of finite size. Depending on
the calculation, the qgrid can be adjusted by specifying it with an
additional flag q_mesh
, for example
vibes output phonopy phonopy/trajectory.son dos q_mesh 26 26 26
The DOS can then be used to evaluate Eq. (6) and (7), i. e., the thermal properties accessible in the harmonic approximation, including the vibrational free energy \(F^\text{ha}\) and the heat capacity at constant volume, \(c_V\). You can compute band structure, density of states, and plot the thermal properties by running:
vibes output phonopy phonopy/trajectory.son full
Carefully inspect all the new files you produced. In the following, you find the DOS plotted next to the band structure (from the file bandstructure_dos.pdf
):
And the figure of the thermal properties (from the file thermal_properties.pdf
)
Animations
At last we highlight that postprocessing files contain animation data that can be visualized. This can be done online, by using amazing phonon website or locally with v_sim
.
In order to view phonon animation online you could follow instructions from FHIvibes Wiki.
If you do full postprocessing of Phonopy output animation files are written in .ascii
format in a subfolder of the output folder. A short introduction on how to use v_sim
to visualize the files is given in the appendix of this tutorial.

https://atztogo.github.io/phonopy/phonopymodule.html?highlight=supercell#supercellmatrix ↩

Note that your results can be a bit different, but the difference should be very small and it is related to numerical accuracy we utilize in these calculations. Also some frequencies could be slightly negative (imaginary), which is also related to the accuracy of the calculations. ↩

W. Setyawan, and S. Curtarolo, Comp. Mat. Sci 49, 299 (2010). ↩↩↩

P. E. Bloechl, O. Jepsen, and O. K. Andersen, Phys. Rev. B. 49, 16223 (1994). ↩

[M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, Group Theory: Application to the Physics of Condensed Matter,Springer (2008).] ↩

L. Chaput, A. Togo, I. Tanaka, and G. Hug, Phys. Rev. B, 84, 094302 (2011). ↩

K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997). ↩

Atsushi Togo and Isao Tanaka, Scr. Mater., 108, 15 (2015). ↩