Phonons from finite differences

First of all to run a phonon calculation a sufficiently large super cell is needed.

The phonon calculations are simply carried out by setting IBRION=5 or IBRION=6 in the INCAR file. IBRION=5, is available as of VASP.4.5, IBRION=6 starting from VASP.5.1. Both flags allow to determine the Hessian matrix (matrix of the second derivatives of the energy with respect to the atomic positions) and the vibrational frequencies of a system.

Important: Only zone centered (Γ-point) frequencies are calculated.

To calculate the Hessian matrix, finite differences are used, i.e. each ion is displaced in the direction of each Cartesian coordinate, and from the forces the Hessian matrix is determined. The two modes differ in the way symmetry is considered. For IBRION=5, all atoms are displaced in all three Cartesian directions, resulting in a significant computational effort even for moderately sized high symmetry systems. For IBRION=6, only symmetry inequivalent displacements are considered, and the remainder of the Hessian matrix is filled using symmetry considerations.

Three parameters influence the determination of the Hessian matrix: The parameter NFREE determines how many displacements are used for each direction and ion, and 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.

Output

The main output ist 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 inc 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 unit cell. The last three nodes are the translational modes (they are usually imaginary modes).

---