Phonon dispersion of bulk Si

Have you set up the local environment?

If not, do that now before proceeding.

The goal of this exercise is to compute the phonon band structure of bulk Si in the diamond structure.

Note

Before computing the vibrational structure, it is always necessary to relax the structure of the system under study, so that the reference configuration truly corresponds to a minimum of the energy, with zero forces.

For this exercise we will just use a lattice constant obtained from a previous optimization with the same basis set and parameters (the value is slightly larger than the experimental one).

Building the supercell to compute the force-constant matrix in real space

Hint

Enter the Si-bulk directory.

If we know the force constant matrix in real space, we can compute the dynamical matrix for every q-point in reciprocal space. For this, we need to create a supercell, in order to get the forces felt by nearby atoms when we displace one of the atoms in the unit cell.

We can use the fcbuild program to generate the supercell that will be used to compute the interatomic force constant matrices in real space. From the time being, we will generate a supercell replicating the unit cell three times (named -1, 0, and 1) along each cartesian direction. In the Si.fcbuild.fdf file you will find these lines:

#
# Options to generate the supercell
#

SuperCell_1    1     # number of shells in which the unit cell is
#   repeated in the direction of the first lattice vector.
SuperCell_2    1     # Idem for the second lattice vector.
SuperCell_3    1     # Idem for the third  lattice vector.

The SuperCell in SIESTA defines the number of shells in which the unit cell is repeated in the direction of the lattice vector. The default value is 0, meaning that the supercell is actually the unit cell.

To generate the supercell, run fcbuild:

fcbuild < Si.fcbuild.fdf

This code dumps the information of the Supercell in an output file, called FC.fdf, that contains

  • The structural data of the supercell, including: number of atoms, lattice constant, lattice vectors, and atomic coordinates and species for all atoms.

  • Additional thing required to compute the interatomic force constants in real space: the indices of the atoms that will be displaced (those belonging to the home unit cell), and the magnitude by which the atoms will be displaced in each of the cartesian directions.

For a final analysis of the results, the supercell should contain enough atoms so that all non-negligible elements of the force constant matrix are computed. The range in real space in which the force constant matrix decays to zero can vary widely from system to system.

Computing the force-constant matrix in real space

Now that we have the supercell info in the FC.fdf file, you can have a look at the Si.ifc.fdf file, which will be the actual input for SIESTA. You can see that some options have “< FC.fdf” next to them, meaning that they will be taken from the FC.fdf file.

The only thing that remains now is to do the actual SIESTA run:

siesta < Si.ifc.fdf > Si.ifc.111.out

The interatomic force constant matrix in real space will be stored in a file called SystemLabel.FC

Computing the q-dependent dynamical matrix, and phonon modes

Once the interatomic force constants in real space have been computed, we need to perfrom a discrete Fourier transform to compute the dynamical matrix in the reciprocal space. Then, the dynamical matrix is diagonalized and its eigenvalues and eigenvectors are computed. This is done using another utility called vibra.

The k-points are defined in the same way as to compute the electronic band structure. We have included these in the same file used to define the supercell, Si.fcbuild.fdf.:

BandLinesScale       pi/a
%block BandLines
1   0.0   0.0   0.0   \Gamma        # Begin at Gamma
45  2.0   0.0   0.0   X             # 45 points from Gamma to X
17  2.0   0.5   0.5   K             # 17 points from X to K
48  2.0   2.0   2.0   \Gamma        # 48 points from K to Gamma
41  1.0   1.0   1.0   L             # 41 points from Gamma to L
%endblock BandLines

Now, you may also notice that this file has additional info in the AtomicCoordinatesAndAtomicSpecies block, namely the atomic mass:

%block AtomicCoordinatesAndAtomicSpecies
-0.125  -0.125  -0.125     1    28.086
 0.125   0.125   0.125     1    28.086
%endblock AtomicCoordinatesAndAtomicSpecies

The atomic mass is necessary for vibra to compute the correct phonon frecuencies. You can then run vibra as:

vibra < Si.fcbuild.fdf

The output of this code includes the following files:

  • SystemLabel.bands: with the different mode frequencies (in cm^-1).

    They are stored in the same way as the electronic band structure.

  • SystemLabel.vectors: with the eigenmodes for each k-points.

Plotting the Phonon Bandstructure

To plot the phonon band structure, proceed in the same way as to plot the electronic band structure, using for example the gnubands utility and later gnuplot:

gnubands  < Si.bands > Si.phonon-bands.111.dat
gnuplot
gnuplot> plot "Si.phonon-bands.111.dat" using 1:2 with lines

You can produce a postscript (image) file of this figure as follows:

gnuplot> set terminal postscript
gnuplot> set output "Si.phonon-bands.111.ps"
gnuplot> replot

Checks of the convergence with the supercell size

One should always check the convergence of the computed phonon band structure with respect the size of the supercell, to be sure that all the relevant interatomic force constant matrix elements are included.

In order to do this:

1. We save all the input and output files we used previously, so that we prevent them from being overwritten:

$ cp Si.fcbuild.fdf Si.fcbuild.111.fdf
$ mv FC.fdf FC.111.fdf
$ mv Si.FC Si.111.FC
$ mv Si.vectors Si.111.vectors
$ mv Si.bands Si.111.bands

2. Edit the file Si.fcbuild.fdf and increase the size of the supercell, adding up to 5 periodic repetitions of the unit cell in each direction (named -2, -1, 0, 1, 2)

SuperCell_1    2     # number of shells in which the unit cell is
#   repeated in the direction of the first lattice vector.
SuperCell_2    2     # Idem for the second lattice vector.
SuperCell_3    2     # Idem for the third  lattice vector.

Note

It might take too long to generate the FC files for the 222 and 333 cases. We provide them in the FILES directory. You will need to adapt the instructions below to this case (rename or copy files as needed).

We then run SIESTA as before, and store the data. For example, for the 2x2x2 case:

fcbuild < Si.fcbuild.fdf
siesta < Si.ifc.fdf > Si.ifc.222.out
vibra < Si.fcbuild.fdf
gnubands < Si.bands > Si.phonon-bands.222.dat
gnuplot
gnuplot> plot "Si.phonon-bands.222.dat" using 1:2 with lines

$ cp Si.fcbuild.fdf Si.fcbuild.222.fdf
$ mv FC.fdf FC.222.fdf
$ mv Si.FC Si.222.FC
$ mv Si.vectors Si.222.vectors
$ mv Si.bands Si.222.bands

And finally, we can compare the results obtained with the three superlattices:

$ gnuplot
gnuplot> plot "Si.phonon-bands.111.dat" using 1:2 with lines,
               "Si.phonon-bands.222.dat" using 1:2 w l,
               "Si.phonon-bands.333.dat" u 1:2 with lines

How well converged was our initial cell?

Note

In case you want some additional background, there is a related slide presentation by Andrei Postnikov.