Quantum-ESPRESSO tutorial: Phonon Dispersions with DFPT

Phonon dispersions with Density-Functional Perturbation Theory

These examples illustrate the applications of Density-Functional Perturbation Theory (DFPT) to the the calculation of phonon dispersions and Interatomic Force Constants in simple semiconductors (Silicon).

Getting ready

A short reminder of the relevant theory can be found here.

The following "make" targets must have been executed:

     make ph pp tools

("make all" is also ok).

You may download the entire exercise file examples_disp.tgz in a directory of your choice. Uncompress and unpack the file and enter in the resulting directory:

     tar -zxvf examples_disp.tgz
     cd examples_disp

In the following, $espresso_dir stands for the root directory of the Quantum-ESPRESSO distribution.

Interatomic Force Constants in Real Space

Whenever the entire phonon dispersion is needed, it is convenient to calculate Interatomic Force Constants (IFC) in real space. This calculation involves three steps:

the usual scf step for the unperturbed system
the calculation of phonons and dynamical matrices on a grid of q-vectors (this can be done in a single step)
the transformation of the dynamical matrices from G- to R-space

Since the range of IFC's in real space is short for nonpolar materials, only a relatively coarse grid of q-vectors in the Irreducible Brillouin Zone is needed. For polar materials, a suitable dipolar term is subtracted out from the IFC's; the remaining analytic term is short-ranged and undergoes the same procedure as in the nonpolar case; finally the non-analytic contribution is added back to the dynamical matrices when needed. Once the IFC's are calculated, it is possible to calculate phonon frequencies and eigenvalues at any wavevector, by just reconstructing the dynamical matrix. If the grid of q-vectors is dense enough, the "reconstructed" dynamical matrix is very close to what one gets by direct calculation, also for a q-vector that was not in the grid of q used to calculate the IFC's in real space. Let us call this a "Fourier interpolation". Coming back to our three steps:

scf step as usual:

     $espresso_dir/bin/pw.x < si.scf.in > si.scf.out

Calculation of dynamical matrices on a grid of q-vectors (will require several minutes).
Sample input file: si.ph.in. One has to specify the grid of q-vectors as in a Monkhorst-Pack grid, with variables "nq1", "nq2", "nq3". The grid must contain Gamma and it must be a subdivision of the reciprocal lattice (i.e. it must correspond to a supercell in R-space). It is automatically generated. One has also to specify "ldisp=.true." to instruct the code to loop over all q-vectors.
```
     $espresso_dir/bin/ph.x < si.ph.in > si.ph.out
```
Note that for each q-vector the dynamical matrices are saved with a different name ("fildyn" + 1-8), while "fildyn" + 0 contains the information on the q-vector grid (type of grid and number of points)
Calculation of IFC's in real space.
This task is performed by auxiliary program "q2r.x". All dynamical matrices are read and Fourier-transformed. The input file for q2r.x contains very little data: file name "fildyn" of dynamical matrices (the same given in input to the phonon code), whether to apply ASR (better to do it: asr='simple' or 'crystal' are ok for simple crystals), the name ("flfrc") of the output file containing the force constants.
```
    $espresso_dir/bin/q2r.x
    &input fildyn='si.dyn', zasr='simple', flfrc='si444.fc' /
```

Phonon Dispersions from Interatomic Force Constants

Once the IFC's in real space are available, one can calculate phonons at any value of q using code "matdyn.x" (do not confuse with "dynmat.x"). It can calculate phonon Density Of States (DOS) as well. Before starting more serious calculations, it is wise to check whether the quality of the Fourier interpolation performed in the previous step is sufficiently good, i.e. if you need a denser q-vector grid (i.e. a longer range for the IFC's) or not. This can be easily checked by comparing the results of matdyn.x and of ph.x for a q-vector that was not in the grid used in the Fourier interpolation. Once you are convinced that your IFC set is ok, let us proceed to two different types of calculations;

A plot of the phonon dispersions along selected high-symmetry lines.
Sample file: si.matdyn.in. You must specify the name of the file ("flfrc") containing the IFC's, and a list of q-points for which the frequencies are to be calculated. You may specify various forms of ASR ('simple' or 'crystal' are ok in this case), where to write the frequencies ("flfrq").
```
     $espresso_dir/bin/matdyn.x < si.matdyn.in
```
The file selected in flfrc, "si.freq", contains a list of frequencies in a format that can be further processed by another auxiliary code, "plotband.x", the same used for band structure plotting:
```
     $espresso_dir/bin/plotband.x si.freq
```
Answer something sensible for Emin and Emax (0 530 for instance; units are cm-1) and for the file in "xmgr" format, for instance, "siph.xmgr"; you may also produce a file directly in printable (postscript) format. The file "siph.xmgr" can be easily plotted using xmgrace or gnuplot plotting programs. Does it look like you expected it?

A plot of the phonon DOS
If you choose option "dos=.true." you must specify the grid of q-vectors used to calculate the DOS (variables "nk1", "nk2, "nk3", using the Monkhorst-Pack logic). You should also specify where you want the DOS written ("fldos") Sample file: si.phdos.in.
```
     $espresso_dir/bin/matdyn.x < si.phdos.in > si.phdos.out
```
The file selected in fldos, "si.phdos", contains the DOS in a format that can be plotted using xmgr, xmgrace, or gnuplot. Can you make the connection between the phonon dispersions with the phonon DOS?

Polar Materials

The calculation of IFC's in polar materials is basically the same as for nonpolar materials. You just have to remember to specify the appropriate option ("epsil=.true.") in the input of trhe phonon code. Can you produce a phonon dispersion plot for AlAs?

Phonons with Ultrasoft pseudopotentials

Phonon calculations with Ultrasoft pseudopotentials are basically the same as for norm-conserving ones. You just have to be aware of the following tricks:

The cutoff for the density, "ecutrho", in the self-consistent calculation, often needs to be quite high (i.e. higher than what is good for total-energy or structural-optimization calculations). In the example provided for Diamond: c.scf.in, a cutoff of 450 Ry is used, while the cutoff for the wavefunction is just 27 Ry.
The threshold for convergence in phonon calculations ("tr2_ph") should be set to quite low values, somewhat lower than what is good for norm-conserving pseudopotentials. In the example provided for the Gamma-point calculation in Diamond: c.ph.in, tr2_ph is set to 1.0E-16.

Can you calculate phonons at Gamma for Diamond?

Phonons in magnetic systems

The example for Ni shows how to calculate phonon dispersions in a magnetic system (incidentally one that is also metallic and with ultrasoft pseudopotentials). The procedure is the same: self-consistency (pw.x) first (sample data: ni.scf.in), followed by the phonon calculation (ph.x) (sample data: ni.ph.in) and by inverse Fourier transform (q2r.x) to get the IFC's in real space (sample data: ni.q2r.in). Once you have the IFC's, you can plot the phonon dispersions using "matdyn.x" (sample data: ni.matdyn.in) and "plotband.x".