Computing the phonon dispersion and DOS

After computing the force constants using the finite differences or density functional perturbation theory approaches, it is possible to compute the phonon dispersion relation as well as the phonon density of states. This is accomplished by Fourier interpolating the interatomic force constants from a supercell calculation to the primitive cell.

Phonon dispersion: Step-by-step instructions

Step 1: Compute the force constants

There are two possible approaches for computing the force constants and then building the dynamical matrix:

1. Using finite differences with (IBRION = 5, 6).
2. Using density functional perturbation theory with (IBRION = 7, 8).

These calculations have to be performed in a supercell so that the force constants vanish at large distances.

Important: If the supercell is not large enough, the phonon frequencies will likely not be properly converged. Always perform convergence studies with respect to the superecell size.

In systems where the unit cell is already sufficiently large, one may directly use the unit-cell geometry.

Step 2: Provide q-points along a high-symmetry path

Create a QPOINTS file containing a q-points path at which the phonon dispersion will be computed. This is accomplished using the line mode of the KPOINTS file format.

Step 3: Compute the phonon dispersion

To plot the phonon dispersion, the tag LPHON_DISPERSION = true must be set in the INCAR file. The amount of information written to the OUTCAR file can be tuned using the (PHON_NWRITE tag).

Reading of force constants

Steps 1-3 can be performed in one VASP calculation. However, generating the finite displacements in the supercell is a time-consuming process. In some situations, it can therefore be beneficial to skip this step and read the force constants from a previous run. To do this, you set LPHON_READ_FORCE_CONSTANTS = true and rename the vaspout.h5 output file from the previous calculation to vaspin.h5. Together with LPHON_DISPERSION = true, this shortcuts the VASP calculation to perform a phonon calculation without having to explicitly calculate the force constants again.

Phonon density of states: Step-by-step instructions

Step 1: Compute the force constants

Same as above. As explained earlier, this can also be skipped if done previously by using LPHON_READ_FORCE_CONSTANTS.

Step 2: Specify a uniform q-point mesh

Create a QPOINTS file that specifies a sufficiently dense, uniform q-point mesh.

Step 3: Compute the density of states

Set PHON_DOS > 0 in the INCAR file. The density of states is computed between $[\omega_{\text{min}}-5\sigma,\omega_{\text{max}}+5\sigma]$ with $\omega_{\text{min}}$ and $\omega_{\text{max}}$ the lowest and highest phonon frequency and $\sigma$ the broadening (PHON_SIGMA).

The number of energy points in this energy range is specified by the PHON_NEDOS tag. To use a Gaussian-smearing method for the computation of the DOS set PHON_DOS = 1 to use the tetrahedron method set PHON_DOS = 2.

Polar materials

If the material is polar, i.e., two or more atoms in the unit cell carry non-zero Born effective charge tensors, the long-range dipole-dipole interaction has to be treated by Ewald summation. This is achieved by setting LPHON_POLAR=.TRUE., supplying the static dielectric tensor (PHON_DIELECTRIC) and the Born-effective charges (PHON_BORN_CHARGES). The values for these have to be obtained from a separate VASP calculation in the unit cell setting the LEPSILON or LCALCEPS INCAR tags. Additionally, the user might specify a reciprocal space cutoff radius (PHON_G_CUTOFF) for the Ewald summation.

Related tags and articles

QPOINTS, LPHON_DISPERSION, PHON_NWRITE, LPHON_POLAR, PHON_DIELECTRIC, PHON_BORN_CHARGES, PHON_G_CUTOFF