Restarting finite differences calculations

From VASP Wiki
Revision as of 13:42, 21 January 2026 by Csheldon (talk | contribs) (→‎Output)
Mind: Available as of VASP 6.6.0

Since VASP 6.6.0, it has been possible to restart finite difference calculations using IBRION = 6 and CHECKPOINT. The displacements are written to a vaspcheck.h5 file. For details of a general finite difference calculation, see the phonons from finite differences HowTo. In this HowTo, we will concern ourselves with restarting and splitting finite difference calculations.

Important: This can only be done using IBRION = 6 currently. We recommend using this generally over IBRION = 5.

Input

There are several options for the CHECKPOINT tag:

  • CHECKPOINT = LEGACY - finite difference calculation are performed as those in VASP 6.5.1. No vaspcheck.h5 is written.
  • CHECKPOINT = RESET (default) - The vaspcheck.h5 file is overwritten or created, and the displacements are written to it.
  • CHECKPOINT = CONTINUE - used for restarting a finite differences calculation from a vaspcheck.h5 file.
  • CHECKPOINT = PREPARE - used for splitting a finite differences calculation.
  • CHECKPOINT = SINGLE - used for running a single displacement after splitting with CHECKPOINT = PREPARE.

We will describe the restart procedure and splitting a calculation below. As an example, we take a 3x3x1 graphene supercell POSCAR file from the phonon tutorials: 

C18
1.0
   7.3521657209830806    0.0000000000000000    0.0000000000000000
  -3.6760828604915403    6.3671622872044793    0.0000000000000000
   0.0000000000000000    0.0000000000000000    8.0000000000000000
C
18
direct
   0.1111111111111133    0.2222222222222200    0.0000000000000000 C
   0.1111111111111133    0.5555555555555532    0.0000000000000000 C
   0.1111111111111133    0.8888888888888866    0.0000000000000000 C
   0.4444444444444466    0.2222222222222200    0.0000000000000000 C
   0.4444444444444466    0.5555555555555532    0.0000000000000000 C
   0.4444444444444466    0.8888888888888866    0.0000000000000000 C
   0.7777777777777801    0.2222222222222200    0.0000000000000000 C
   0.7777777777777799    0.5555555555555532    0.0000000000000000 C
   0.7777777777777799    0.8888888888888866    0.0000000000000000 C
   0.2222222222222200    0.1111111111111133    0.0000000000000000 C
   0.2222222222222200    0.4444444444444466    0.0000000000000000 C
   0.2222222222222199    0.7777777777777799    0.0000000000000000 C
   0.5555555555555532    0.1111111111111133    0.0000000000000000 C
   0.5555555555555534    0.4444444444444466    0.0000000000000000 C
   0.5555555555555532    0.7777777777777799    0.0000000000000000 C
   0.8888888888888866    0.1111111111111133    0.0000000000000000 C
   0.8888888888888866    0.4444444444444466    0.0000000000000000 C
   0.8888888888888866    0.7777777777777799    0.0000000000000000 C

along with a 4x4x1 k-mesh in our KPOINTS file: 

K points
 0
Gamma
4  4  1
0  0  0

and PAW C_s 04May1998 POTCAR. We can use the following INCAR file: 

SYSTEM = graphene
ENCUT = 400

# electronic
PREC = Accurate
NELMIN = 5
EDIFF = 1e-8
ISMEAR = -1
SIGMA = 0.2
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.

# ionic (finite differences)
IBRION = 6
POTIM  = 0.015

Output

We assume that this calculation has failed after a certain point or accidentally been cancelled: 

DAV:   9    -0.181350430125E+03   -0.25600E-06   -0.38452E-08   848   0.118E-03    0.774E-04
DAV:  10    -0.181350430581E+03   -0.45554E-06   -0.20686E-08   864   0.738E-04    0.208E-04
DAV:  11    -0.181350430676E+03   -0.94857E-07   -0.22346E-09   704   0.354E-04    0.207E-04
srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
slurmstepd-test01: error: *** JOB 254054 ON test01 CANCELLED AT 2026-01-19T16:25:49 ***
slurmstepd-test01: error: *** STEP 254054.0 ON test01 CANCELLED AT 2026-01-19T16:25:49 ***

If you look in the directory, you will see the vaspcheck.h5 file. Inside it will contain the displacement calculations that have been completed up to the point of the crash: 

h5ls vaspcheck.h5
data-1                   Group
displacements            Group
symmetry                 Group

The calculation can then be restarted by adding CHECKPOINT = CONTINUE to the INCAR:


The phonon modes and frequencies are written to the OUTCAR file after 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 Cartesian coordinates.

There should be 3[math]\displaystyle{ N }[/math] normal modes, where [math]\displaystyle{ N }[/math] 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).

Finally, IBRION=6 and ISIF≥3 allows to calculate the elastic constants. The elastic tensor is determined by performing six finite distortions of the lattice and deriving the elastic constants from the strain-stress relationship.[1] The elastic tensor is calculated both, for 'clamped' ions, as well, as allowing for relaxation of the ions. The elastic moduli for rigid ions are written after the line 

SYMMETRIZED ELASTIC MODULI (kBar)

The ionic contributions are determined by inverting the ionic Hessian matrix and multiplying with the internal strain tensor,[2] and the corresponding contributions are written after the lines: 

ELASTIC MODULI CONTR FROM IONIC RELAXATION (kBar)

The final elastic moduli, including both, the contributions for distortions with rigid ions and the contributions from the ionic relaxations, are summarized at the very end: 

TOTAL ELASTIC MODULI (kBar)

There are a few caveats to this approach: most notably, the plane-wave cutoff (ENCUT) needs to be sufficiently large to converge the stress tensor. This is usually only achieved if the default cutoff is increased by roughly 30%, but it is strongly recommended to increase the cutoff systematically, (e.g., in steps of 15%), until full convergence is achieved.

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). Additionally, compare the result to phonons from density-functional-perturbation theory (DFPT).

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

grep THz OUTCAR

To get an accurate phonon dispersion, perform the force-constants calculation in a large enough supercell. When increasing the size of the supercell, we recommend decreasing the k-point density in the KPOINTS file to yield the same resolution. 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[3] to post-process the results of a finite differences calculation done with VASP.[4]

Tip: In contrast to computing phonons within DFPT, the finite difference approach can be used in combination with any Exchange-correlation functional.

IBRION=5, is available as of VASP.4.5, IBRION=6 starting from VASP.5.1. In some older versions (pre VASP.5.1), NSW (number of ionic steps) must be set to 1 in the INCAR file, since NSW=0 sets the IBRION=−1 regardless of the value supplied in the INCAR file. Although VASP.4.6 supports IBRION=5-6, VASP.4.6 does not change the set of k points automatically (often the displaced configurations require a different k-point grid). Hence, the finite difference routine might yield incorrect results in VASP.4.6.

Related tags and sections

IBRION, ISIF, POTIM,

Phonons: Theory

Phonons from density-functional-perturbation theory, Computing the phonon dispersion and DOS

References

  1. ↑ Y. Le Page and P. Saxe, Phys. Rev. B 65, 104104 (2002).
  2. ↑ X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B 72, 035105 (2005).
  3. ↑ http://phonopy.github.io/phonopy/index.html (2022).
  4. ↑ http://phonopy.github.io/phonopy/vasp-dfpt.html (2022).
Retrieved from "https://vasp.at/wiki/index.php?title=Restarting_finite_differences_calculations&oldid=33320"