Exercise 2: Running Phonopy calculations via FHI-vibes

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 FHI-vibes 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:                       1e-05
    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 FHI-vibes 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 ¹:

\[ \begin{pmatrix} \boldsymbol{a}_\text{s} & \boldsymbol{b}_\text{s} & \boldsymbol{c}_\text{s} \end{pmatrix} = \begin{pmatrix} \boldsymbol{a}_\text{u} & \boldsymbol{b}_\text{u} & \boldsymbol{c}_\text{u} \end{pmatrix} \cdot M_\text{S}~. \label{eq:smatrix} \]

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 ⁸, force constants are generated based on the finite displacement method ⁶. 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 Parlinski-Li-Kawazoe method ⁷.

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, FHI-vibes 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\))²! 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 FHI-vibes 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. ³

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) ⁵.

You can list the high symmetry points of the lattice of your geometry with FHI-vibes 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 so-called 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 ³, as shown above.

In order to compute the phonon band structure, we perform the following steps:

• Run the Phonopy calculations with FHI-vibes 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 post-processing 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 post-processing, 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 so-called Tetrahedron method to interpolate between the points ⁴. Afterwards it literally counts the number of frequencies in bins of finite size. Depending on the calculation, the q-grid 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 post-processing 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 FHI-vibes 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.

---

1. https://atztogo.github.io/phonopy/phonopy-module.html?highlight=supercell#supercell-matrix ↩

2. 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. ↩

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

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

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

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

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

8. Atsushi Togo and Isao Tanaka, Scr. Mater., 108, 1-5 (2015). ↩