Revision as of 16:27, 4 July 2024

Band-structure calculations find the energy of a cell using a finite number of plane waves and a finite number of k-points. When comparing cells of different sizes, i.e. volume, this results in them each having different plane wave basis sets. This would be solved by using an infinite number of k-points and plane waves. In practice, a large enough plane wave energy cutoff and number of k-points leads to converged energies.^[1] However, when the basis set is too small, i.e. prematurely truncated, this results in discontinuities in the total energy between cells of varying volumes. These discontinuities between energy and volume create stress that decreases the volume compared to fully converged calculations, due to the diagonal components of the stress tensor being incorrect.^[2] This is the "Pulay stress". Figures 1 and 2 show these discontinuities for diamond for prematurely truncated energy cutoffs and k-point meshes, respectively, in comparison to converged curves. It may be corrected by using a few different methods in VASP. It is important to note that problems due to the Pulay stress can often be neglected if only volume conserving relaxations are performed. This is because the Pulay stress is, usually, nearly uniform and only changes the diagonal elements of the stress tensor by a constant amount.

How to correct

The Pulay stress may be corrected in multiple ways. Generally, by calculating the relaxed structure with a larger basis set by increasing the ENCUT until convergence is reached:

Set ENCUT to $1.3\times$ the default cutoff, or PREC=High in VASP.4.4. N.B. this is the lower limit.
Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.
Further increase the ENCUT and repeat Steps 1 and 2, until the structure no longer changes, i.e. is converged.

If volume relaxations are necessary, the following two procedures may be followed:

1. Volume relaxation

One way is to very accurately relax the structure in a series of calculations on one structure:

Step 1.

Relax from the starting structure with Gaussian or Methfessel-Paxton smearing (ISMEAR = 0 or 1).

Step 2.

Copy the CONTCAR to POSCAR and relax the structure again.

Step 3.

Change the smearing method to the tetrahedron method (i.e. ISMEAR=-5) and perform a single point calculation, i.e. no relaxation of structure. These are will give highly accurate energies.

(Optional)

The previous steps should yield good energies. If there are still problems, such as the Energy vs. Volume curve remaining jagged, a few more additional steps may be tried:

Step 4a.

Try further increasing the ENCUT. Alternatively, improve the FFT grid by setting PREC=Accurate.

Step 4b.

To avoid additional computational cost due to increased cutoff energy, the STRESS output in VASP may be corrected by subtracting PSTRESS. The Pulay stress is only weakly dependent on volume and ionic configuration; it is mainly determined by the composition. In addition, all volume relaxations will take PSTRESS into account. It is important keep in mind that PSTRESS should only be used if increasing the cutoff is not a viable option.

Perform a single point calculation for starting structure with a higher ENCUT. A good estimation for the Pulay stress is given as a negative pressure in the OUTCAR file:

  external pressure =    -100.29567 kB

Add PSTRESS to the INCAR.
Set it equal to the Pulay stress given by the 'external pressure' in the OUTCAR.

2. Equation of state fitting

An alternative way to avoid relaxing the volume is to relax the ionic positions and cell shape for a fixed set of volumes, i.e. multiple POSCAR files. These are then fitted to an equation of state, e.g. Murnaghan.^[3] As the Pulay stress is almost isotropic, only a constant value is added to the diagonal elements of the stress tensor. Therefore, the relaxation for a fixed volume will yield highly accurate structures. This approach has the advantage of also providing the bulk modulus, and we have found it may be used safely with the default energy cutoff.

The procedure is similar to the previous section. In this case, one has to do the calculations for a set of fixed volumes. The following steps have to be done in these calculations:

Step 1.

Select one cell and relax it from the starting structure, keeping the volume fixed while the ionic positions and the cell shape are free to change (ISIF=4; ISMEAR=0 or 1).

Step 2.

Copy the CONTCAR to POSCAR and relax the structure again (if the initial cell shape was reasonable this step can be skipped if the cell shape is kept fixed, you never have run VASP twice).

Step 3.

Change the smearing method to the tetrahedron method (i.e. ISMEAR=-5) and perform a single point calculation, i.e. no relaxation of structure. These will give highly accurate energies.

Step 4.

Repeat Steps 1-3 for as many different volumes as desired and then fit the Energy vs. Volume curve to your chosen Equation of State.

Possible issues and advice on how to address them

When producing Energy vs. Volume plots, improper settings may result in jagged curves, cf. Fig. 1 and 2. This is commonly due to two reasons:

Basis set incompleteness

The basis set is discrete and incomplete, so when the volume changes, additional plane waves are added. That causes small discontinuous changes in the energy.
1. Use a larger plane wave cutoff, cf. Fig. 1. This is usually the preferred and simplest solution.
2. Use more k-points, cf. Fig. 2.^[4]
Discontinuity of FFT grids

Different precisions of the FFT grid defined by PREC may be used, e.g. Normal or Accurate.^[5] This introduces small errors when the volume changes, as the FFT grids change discontinuously. In other words, at each volume a different FFT grid is used, causing the energy to jump discontinuously between different volumes. For more details on FFT grids, see ^[6].
1. Use PREC=Accurate
2. Increase the plane wave cutoff.
3. Set your FFT grids manually, and choose the one that is used per default for the largest volume.

References

↑ G. P. Francis and M. C. Payne, Finite basis set corrections to total energy pseudopotential calculations, J. Condens. Matter Phys. 2, 4395 (1990).
↑ All volume and cell shape relaxation algorithms in VASP work with a constant basis set. In this way, all energy changes are strictly consistent with the calculated stress tensor. This results in an underestimation of the equilibrium volume, unless a large plane wave cutoff is used. By keeping the basis set constant during relaxations, all the initial G-vectors within a sphere are included in the basis. As the cell shape is relaxed, the direct and reciprocal lattice vectors change. The number of reciprocal G-vectors in the basis is kept fixed but the length of the G-vectors changes, indirectly changing the energy cutoff. In other words, the shape of the cutoff region changes from a sphere to an ellipsoid. Restarting VASP after a volume relaxation causes a new "spherical" cutoff to be adopted, thus changing the energy discontinuously.
↑ Murnaghan Equation of State, www.wikipedia.org (2024)
↑ This solves the problem because the criterion for including a plane wave in the basis set is $\vert {\bf G} + {\bf k} \vert < {\bf G}_{\rm cut}$ . This means at each k-point a different basis set is used, and additional plane waves are added at each k-point at different volumes. In turn, the energy vs. volume curve becomes smoother.
↑ For PREC=Accurate, the FFT grids are chosen such that ${\bf H} \vert \phi>$ is exactly evaluated. Whereas, for PREC=Normal the FFT grids are set to 3/4 of the value that is required for an exact evaluation of ${\bf H} \vert \phi>$ .
↑ Introduction to ab-initio simulation in VASP

[payne:francis:1990-1] G. P. Francis and M. C. Payne, Finite basis set corrections to total energy pseudopotential calculations, J. Condens. Matter Phys. 2, 4395 (1990).

[2] All volume and cell shape relaxation algorithms in VASP work with a constant basis set. In this way, all energy changes are strictly consistent with the calculated stress tensor. This results in an underestimation of the equilibrium volume, unless a large plane wave cutoff is used. By keeping the basis set constant during relaxations, all the initial G-vectors within a sphere are included in the basis. As the cell shape is relaxed, the direct and reciprocal lattice vectors change. The number of reciprocal G-vectors in the basis is kept fixed but the length of the G-vectors changes, indirectly changing the energy cutoff. In other words, the shape of the cutoff region changes from a sphere to an ellipsoid. Restarting VASP after a volume relaxation causes a new "spherical" cutoff to be adopted, thus changing the energy discontinuously.

[murnaghan:web-3] Murnaghan Equation of State, www.wikipedia.org (2024)

[4] This solves the problem because the criterion for including a plane wave in the basis set is $\vert {\bf G} + {\bf k} \vert < {\bf G}_{\rm cut}$ . This means at each k-point a different basis set is used, and additional plane waves are added at each k-point at different volumes. In turn, the energy vs. volume curve becomes smoother.

[5] For PREC=Accurate, the FFT grids are chosen such that ${\bf H} \vert \phi>$ is exactly evaluated. Whereas, for PREC=Normal the FFT grids are set to 3/4 of the value that is required for an exact evaluation of ${\bf H} \vert \phi>$ .

[vasp:intro:lecture:web-6] Introduction to ab-initio simulation in VASP

[1]

[2]

[3]

[4]

[5]

[6]

@@ Line 30: / Line 30: @@
 === Step 4b. ===
-To avoid additional computational cost due to increased cutoff energy, the ''STRESS'' output in VASP may be corrected by subtracting {{TAG|PSTRESS}}. In addition, all volume relaxations will take {{TAG|PSTRESS}} into account. It is important keep in mind that {{TAG|PSTRESS}} should only be used if increasing the cutoff is not a viable option.
+To avoid additional computational cost due to increased cutoff energy, the ''STRESS'' output in VASP may be corrected by subtracting {{TAG|PSTRESS}}. The Pulay stress is only weakly dependent on volume and ionic configuration; it is mainly determined by the composition. In addition, all volume relaxations will take {{TAG|PSTRESS}} into account. It is important keep in mind that {{TAG|PSTRESS}} should only be used if increasing the cutoff is not a viable option.
 <ol type = 'i'>
 <li>Perform a single point calculation for starting structure with a higher {{TAG|ENCUT}}. A good estimation for the Pulay stress is given as a negative pressure in the {{TAG|OUTCAR}} file:</li>