Phonons from finite differences

The phonon calculations using a finite differences approach are carried out by setting IBRION=5 or IBRION=6 in the INCAR file. These flags allow determining the second-order force constants matrices and the Hessian matrix (matrix of the second derivatives of the energy with respect to atomic displacements) and the vibrational frequencies of a system. If ISIF>=3 the internal strain tensors are computed as well.

Mind: Only zone-center (Γ-point) frequencies are calculated.

It is possible to obtain the phonon dispersion at different q points by computing the force constants on a sufficiently large supercell and Fourier interpolating the dynamical matrices in the primitive cell.

Theory

The second-order force constants are computed using finite differences of the forces when each ion is displaced in each independent direction. This is done by creating systems with finite ionic displacement of atom $a$ in direction $i$ with magnitude $\lambda$, computing the orbitals $\psi _{\lambda }^{u_{i}^{a}}$ and the forces for these systems. The second-order force constants are then computed using

$\Phi _{ij}^{ab}={\frac {\partial ^{2}E}{\partial u_{i}^{a}\partial u_{j}^{b}}}=-{\frac {\partial F_{i}^{a}}{\partial u_{j}^{b}}}\approx -{\frac {\left(\mathbf {F} [\{\psi _{\lambda /2}^{u_{j}^{b}}\}]-\mathbf {F} [\{\psi _{-\lambda /2}^{u_{j}^{b}}\}]\right)_{i}^{a}}{\lambda }},\quad {a=1,..,N_{\text{atoms}}}\quad {b=1,..,N_{\text{atoms}}}\quad {i=x,y,z}\quad {j=x,y,z}$

where $u_{i}^{a}$ corresponds to the displacement of atom $a$ in the cartesian direction $i$ and $\mathbf {F} [\psi ]$ retrieves the set of forces acting on all the ions given the $\psi _{n\mathbf {k} }$ orbitals.

Similarly, the internal strain tensor is

$\Xi _{il}^{a}={\frac {\partial ^{2}E}{\partial u_{i}^{a}\partial u_{j}^{b}}}={\frac {\partial \sigma _{l}}{\partial u_{i}^{a}}}\approx {\frac {\left(\mathbf {\sigma } [\{\psi _{\lambda /2}^{u_{i}^{a}}\}]-\mathbf {\sigma } [\{\psi _{-\lambda /2}^{u_{i}^{a}}\}]\right)_{l}}{\lambda }},\qquad {l=xx,yy,zz,xy,yz,zx}\quad {i=x,y,z}\quad {a=1,..,N_{\text{atoms}}}$

where $\mathbf {\sigma } [\psi _{n\mathbf {k} }]$ computes the strain tensor given the $\psi _{n\mathbf {k} }$ orbitals.

Input

There are two options to compute the second order force-constants using finite differences:

• IBRION=5, all atoms are displaced in all three Cartesian directions, resulting in a significant computational effort even for moderately sized high-symmetry systems.
• IBRION=6, only symmetry inequivalent displacements are considered, and the remainder of the force constants matrix is filled using symmetry operations.

POTIM determines the step size. The step size is defaulted to 0.015 Å (starting from VASP.5.1), if too large values are supplied in the input file. Expertise shows that this is a very reasonable compromise.

The NFREE tag determines how many displacements are used for each direction and ion:

• NFREE=2 uses central differences, i.e., each ion is displaced by a small positive and negative displacement, ±POTIM, along each of the Cartesian directions.
• NFREE=4 uses four displacements along each of the Cartesian directions ±POTIM and ±2×POTIM.
• NFREE=1 uses a single displacement (this is strongly discouraged since it can ).

Output

The main output is written to the OUTCAR file and starts with the following lines:

 Eigenvectors and eigenvalues of the dynamical matrix
 ----------------------------------------------------

The following lines are repeated for each normal mode and a should look like the following example output:

   1 f  =   14.329944 THz    90.037693 2PiTHz  477.995462 cm-1    59.263905 meV
             X         Y         Z           dx          dy          dz
      0.000000  0.000000  0.000000     0.009046   -0.082007   -0.006117
      0.000000  2.731250  2.731250     0.009046    0.106244    0.006563
      0.000000  5.462500  5.462500     0.009046    0.082007    0.006117
      0.000000  8.193750  8.193750     0.009046   -0.106244   -0.006563
      ...
   2 f  =   14.329944 THz    90.037693 2PiTHz  477.995462 cm-1    59.263905 meV
             X         Y         Z           dx          dy          dz
      0.000000  0.000000  0.000000     0.003458    0.021825   -0.093181
      0.000000  2.731250  2.731250     0.003458    0.005416    0.094689
      0.000000  5.462500  5.462500     0.003458   -0.021825    0.093181
      0.000000  8.193750  8.193750     0.003458   -0.005416   -0.094689
      ...
   ...

The first number is the number for the normal mode. If it is followed by f, it is a mode on the real axis (vibrationally stable). Otherwise if it is followed by f/i, the mode is an imaginary mode ("soft mode"). The following entries are just the eigenfrequency of the mode in different units.

The following column is the label for the atomic positions in Cartesian coordinates (x,y,z) and the normalized eigenvectors of the eigenmodes in direct coordinates.

There should be 3$N$ normal modes, where $N$ is the number of atoms in the supercell (POSCAR). The modes are ordered in descending order with respect to the eigenfrequency. The last three modes are the translational modes (they are usually disregarded).

Practical hints

To get the phonon frequencies quickly on the command line simply type the following:

grep THz OUTCAR

---