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. When these flags are set the second-order force constants are computed using finite differences. The dynamical matrix is constructed, diagonalized and the phonon modes and frequencies of the system are reported in the OUTCAR file. 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 second-order force constants on a sufficiently large supercell and Fourier interpolating the dynamical matrices in the unit cell.

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. If too large values are supplied in the input file the step size is defaulted to 0.015 Å (starting from VASP.5.1). 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).

Selective dynamics are presently only supported for IBRION=5; in this case, only those components of the Hessian matrix are calculated for which the selective dynamics tags are set to .TRUE. in the POSCAR file. Contrary to the conventional behavior, the selective dynamics tags refer to the Cartesian components of the Hessian matrix. For the following POSCAR file, for instance,

Cubic BN
   3.57
 0.0 0.5 0.5
 0.5 0.0 0.5
 0.5 0.5 0.0
   1 1
selective
Direct
 0.00 0.00 0.00  F F F
 0.25 0.25 0.25  T F F

atom 2 is displaced in the x-direction only, and only the x-component of the second atom of the Hessian matrix is calculated.

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 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 label of the normal mode. If this number is followed by f it is a purely real mode, stating the mode is vibrationally stable. Otherwise, if it is followed by f/i, the mode is an imaginary mode ("soft mode"). These labels are followed by the eigenfrequency of the mode in different units.

The following table labeled by (x,y,z,dx,dy,dz) contains the Cartesian positions of the atoms 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

The computation of the second-order force constants requires accurate forces. Therefore, the tag PREC=Accurate is recommended in the INCAR. The ADDGRID=TRUE should not be set without careful testing.

A practical way to check for convergence is to monitor the Γ point (q=0) optical mode frequencies while varying the ENCUT, PREC, and k point density (KPOINTS). This should be done in the unit cell POSCAR while setting IBRION=6 or 7 in your INCAR.

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

grep THz OUTCAR

To get an accurate phonon dispersion it is important to perform the force-constants calculation in a large enough supercell. When increasing the size of the supercell it is recommended to decrease the k point density in the KPOINTS file. For example, for the primitive cell of silicon, a 12x12x12 Gamma-centered k point mesh is needed to obtain accurate phonon frequencies at the Gamma point. When replicating the unit cell to a 2x2x2 supercell a 6x6x6 k point mesh will produce an equivalent sampling. For a 4x4x4 supercell a 3x3x3 k point mesh will suffice.

It is possible to use phonopy to post-process the results of a finite differences calculation done with VASP.