# Advanced topics in phonons¶

- author
Andrei Postnikov (Univ. Lorraine)

Note

For background, see the excellent slide presentation for a lecture by Andrei Postnikov in a previous Siesta school.

Note

There is some extra material for this tutorial here.

Note

To fully appreciate the issues in this tutorial you need to have understood the concepts presented in the basic vibrations tutorial.

## The supercell explosion¶

We are going to consider a monolayer of BN, with two atoms per unit
cell. In directory *Phonon_Q* we have already generated a 5x5x1
supercell with *fcbuild* (see the basic vibrations tutorial), by
specifying:

```
Supercell_1 2
Supercell_2 2
Supercell_3 0
```

The supercell contains 50 atoms, and we should do in principle 2x6=12 calculations displacing the two atoms in the unit cell in the positive and negative cartesian directions to compute the forces on all the atoms in the supercell. That is a lot of work, so we have prepared a FC file for the rest of the exercise.

You need to copy this pre-generated FC file: `cp BN.FC_save BN.FC`

,
and then run vibra: `vibra < BN_vibra.fdf`

. You will get the *BN.bands*
and *BN.vectors* files.

The *BN.bands* file holds the phonon dispersion information, that you
can plot with:

```
gnubands -G -o bnq BN.bands
gnuplot --persist bnq.gplot
```

You will get something like

which looks quite bad in a neighborhood of Gamma.

Obviously, to improve on this we should increase the size of the supercell, maybe going to 7x7x1, 9x9x1, or even larger. The prospect is not pleasing.

But there is another way. Instead of generating a “brute force” supercell, we can generate a supercell that will describe correctly the dispersion along the high-symmetry directions that we care about. In this particular case, we can construct, following the indications in the slide presentation, a supercell that will exactly map the K and M points to Gamma. This is not the supercell for bookkeeping of the force-constants, but a new “unit cell” that serves our purposes.

The process starts in the BZ, chosing new unit vectors to define a lattice that will include M and K:

The reciprocal cell of this new unit cell is the real-space supercell we seek:

In directory *Phonon_GMK* we have the necessary files to compute the
phonons **at Gamma** in this new “enlarged unit cell”. Now, with 24
atoms, we would need 24*6=144 calculations. It is a bit more CPU than
the original case above (12 calculations for 50 atoms, each about 8
times more expensive than one for 24 atoms). Instead of running the
calculations, you can use the provided FC file. (Copy the *BN_save.*
files to use simply the prefix *BN*)

Run `vibq0`

and use the *BN* SystemLabel when asked. It will read the
BN.FC and BN.XV files and produce
*.bands* and *.vectors* files, just like `vibra`

, but with a much
simpler interface when only the Gamma point is handled.

In this case, the plotting is different because there is a single q-point:

```
gnubands -G -o gamma BN.bands
gnuplot
gnuplot> plot "gamma" using 1:2 with points
```

and you will get a column of points marking the frequencies of the
zone-center phonons in this enlarged unit cell. Not very useful
compared to the original dispersion. All the modes are piled up at
Gamma, and we need to **unfold** them.

Note

There are fundamental similarities between this problem and that of the unfolding of electronic bands covered in this tutorial, but the techniques are different.

### Unfolding the Gamma “phonon structure” in the supercell¶

The character of these different modes can be revealed by projecting their eigenvectors onto waves with different Q-vector (see the lecture liked-to above):

The projection is done by the `phdos`

program which is part of the
lattice dynamics analysis package.

In order to probe projections onto different symmetric points of the BZ, the following coordinates of Q-points can be specified (in Cartesian coordinates and units of 1/Bohr, omitting the 2*pi factor). For A=2.495 Bohr:

```
M [0, 2*pi/(a*sqrt(30), 0] -> [0, 0.2314029, 0], or
M [pi/a, pi/(a*sqrt(3)), 0] -> [0.2004008, 0.1157014, 0];
K [2*pi/(3*a), 2*pi/(a*sqrt(3)), 0] -> [0.1336005, 0.2314029, 0];
K [4*pi/(3*a), 0, 0] -> [0.2672011, 0, 0]
```

according to the picture:

We have prepared the *G.inp*, *K.inp*, *M.inp* files to run the
program. For example:

```
phdos < G.inp
mv BN.WSQ BN.G.WSQ
mv BN.VSQ BN.G.VSQ
```

will produce a *.VSQ* and a *.WSQ* file, which we rename for
convenience and clarity. To plot the contents of the
second one, which is more useful:

```
gnuplot
gnuplot> plot "BN.G.WSQ" u 1:3 w l
gnuplot> plot "BN.G.WSQ" u 1:2 w l
```

We will get sharp peaks at the positions of the “true” Gamma phonons (three acoustic modes with zero frequency, an intermediate one, and a degenerate optical mode), as in this figure:

We can do the same for the K and M points. For K:

For M:

The two different colors correspond to the B and N atoms’ role in each mode. Being lighter, B will feature more prominently in the higher-frequency modes.

Can you check whether the phonons at K and M match the values in the general dispersion in the first section?

Can you think of which q-vectors to use with `phdos`

to see whether
we get better results in the lines Gamma-K and Gamma-M?

Note

One would be tempted to get a “dispersion” by selecting arbitrary
q-vectors along a line and feeding them to `phdos`

. However, it
should be noted that we have really computed only 3N modes, where N
is the number of atoms in the supercell. For some *q*’s we will get
sharp peaks, but for others the signal will degrade and it would be
difficult to make out “frequencies”.

Hint

The *.VSQ* files might be useful to get the frequencies of the peaks,
as they have no broadening.

Note

There is an alternative method to compute the phonon dispersion at
arbitrary *q*s in the BZ, using only *unit cell* calculations. It
is based on Density-Functional Perturbation Theory. While
convenient, this route removes half the fun of understanding
supercells and folding issues.

## LO-TO splitting and Born dynamical charges¶

Here is a picture of the phonon dispersion in hexagonal BN:

The LO-TO splittings at Gamma, in a *3D* system, are caused by the effect of the electric
field that the optical modes produce. There is an extra term in the
force constants:

The second term is non-analytic, in the sense that it has different limiting values when q–>0 depending on the direction in which we approach it.

The *Z*’s in the above equation are the Born effective or dynamical
charges, and can be computed as the derivative of the polarization in
the crystal with respect to the atomic displacement, as shown in
this tutorial, which you might enjoy
following. In fact, the Born charges for hexagonal BN are explicitly
calculated there, so you could go ahead and estimate roughly the size
of the extra terms in the dynamical matrix for specific directions.

There is one extra ingredient: the \(\varepsilon_\infty\) in the non-analytic term is the clamped-ion dielectric constant, which can be computed using the methods outlined in the tutorial on optical properties.

We will not carry out any more specific calculations here, but it is worth noting the subtle links of the topic of phonons with other domains.

Also, it is apparent that to fully analyze a given system you might
have to do several kinds of simulations: force-constants, Born
charges, dielectric constant, etc. Ideally, you would set up an automated
**workflow** that can help you manage everything. See the AiiDA
framework tutorial for more information.

Note

Another use of the Born charges is in computing the infrared intensities of the modes. Only a few modes (those which involve changes in the polarization) would be infrared-active.

The `vibra`

program can compute the infrared intensities if both
the FC and the BC (for Born charges) files are available. See the
Siesta and Vibra manuals for more information.