Pulay stress

Pseudopotential calculations are performed to calculate the energy of a cell using a finite number of plane waves and a finite number of k-points. When comparing between cells of different sizes, 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. However, when the basis set is too small, i.e. prematurely truncated, this results in discontinuities in the total energy between cells of varying size, i.e. varying volume. These discontinuities between energy and volume creates stress that decreases the volume compared to fully converged calculations, due to the diagonal components of the stress tensor being incorrect. This is "Pulay stress".

Pulay Stress

How to calculate the Pulay stress

The Pulay stress shows only a weak dependency on volume and the ionic configuration. It is mainly determined by the composition. The simplest way to estimate the Pulay stress is to relax the structure with a large basis-set ( $1.3\times$ default cutoff is usually sufficient, or PREC=High in VASP.4.4). Then re-run VASP for the final relaxed positions and cell parameters with the default cutoff or the desired cutoff. Look for the line 'external pressure' in the OUTCAR file:

 external pressure =    -100.29567 kB

The corresponding (negative) pressure gives a good estimation of the Pulay stress.

Accurate bulk relaxations with internal parameters (one)

The general message is: whenever possible avoid volume relaxation at the default energy cutoff. Either increase the basis set by setting ENCUT manually in the INCAR file, or use the method two suggested below. This avoids doing volume relaxations at all. If volume relaxations are the only possible and feasible option please use the following step by step procedure (which minimizes errors to a minimum):

Relax from starting structure (ISMEAR should be 0 or 1).
Start a second relaxation from previous CONTCAR file (re-relaxation).
As a final step, perform one more energy calculation using the tetrahedron method switched on (i.e. ISMEAR=-5), to obtain highly accurate energies (no relaxation for the final run). Possibly increase the energy cutoff even further.

A few things should be remarked here: never used the energy obtained at the end of a relaxation run, if you allow for cell shape relaxations (the final basis set might not correspond to the desired spherical cutoff sphere). Instead, perform one additional static run after completing the relaxation.

If the relaxation will yield a structure with reasonably small structural "errors", the final error in the energy of step 3 is only of second-order (with respect to the structural errors). If you take the energy directly from the relaxation run, errors are usually significantly larger. Another important point is that the most reliable results for the relaxation are obtained if the starting cell parameters are very close to the final cell parameters. If different runs yield different results, then the one run that started from the configuration, which was closest to the relaxed structure, is the most reliable one.

We strongly recommend doing any volume (and to a lesser extent cell shape) relaxation with an increased basis set. ENCUT= $1.3\times$ default cutoff is reasonable accurate in most cases. PREC=High does also increase the energy cutoff by a factor of 1.25. At an increased cutoff the Pulay-stress correction is usually not required.

Finally, if the default cutoff is used for the relaxation, the PSTRESS tag should be set in the INCAR file: evaluate the Pulay stress along the guidelines given in the previous section and add an input line to the INCAR file reading (usually a negative number):

PSTRESS = Pulay stress

From now on all STRESS output of VASP is corrected by simply subtracting PSTRESS. In addition, all volume relaxations will take PSTRESS into account. It should be said again, use PSTRESS only if increasing the cutoff is not a viable option for some reason.

Accurate bulk relaxations with internal parameters (two)

It is possible to avoid volume relaxation in many cases: The method we have used quite often in the past is to relax the structure (cell shape and internal parameters) for a set of fixed volumes (ISIF=4). The final equilibrium volume and the ground-state energy can be obtained by a fit to an equation of state. The reason why this method is better than volume relaxation is that the Pulay stress is almost isotropic, and thus adds only a constant value to the diagonal elements of the stress tensor. Therefore, the relaxation for a fixed volume will yield highly accurate structures.

The outline for such a calculation is almost the same as in the previous section. But in this case, one has to do the calculations for a set of fixed volumes. At first sight, this seems to be much more expensive than method number one (outlined in the previous section). But in many cases, the additional costs are only small, because the internal parameters do not change very much from volume to volume. The following steps have to be done in these calculations:

Select one volume and relax from starting structure keeping the volume fixed (ISIF=4; ISMEAR=0 or 1).
Start a second relaxation from the previous CONTCAR file (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).
As a final step, perform one more energy calculation with the tetrahedron method switched on (ISMEAR=-5), to get very accurate unambiguous energies (no relaxation for the final run).

The method has also other advantages, for instance, the bulk modulus is readily available. We have found in the past that this method can be used safely with the default cutoff.

FAQ: Why is my energy vs. volume plot jagged

This is a very common question from people who start to do calculations with plane wave codes. There are two reasons why the energy vs. volume plot looks jagged:

Basis set incompleteness. The basis set is discrete and incomplete, and when the volume changes, additional plane waves are added. That causes small discontinuous changes in the energy.
- Solutions:
  - Use a larger plane wave cutoff: This is usually the preferred and simplest solution.
  - Use more k-points : 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.
A second possible reason for the jagged E(V) curve is the use of PREC=Normal. For PREC=Accurate the FFT grids are chosen such that ${\bf H} \vert \phi>$ ${\bf H} \vert \phi>$ is exactly evaluated. 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>$ ${\bf H} \vert \phi>$ . This introduces small errors because when the volume changes 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.
- Solutions:
  - Use PREC=Accurate, or increase the plane wave cutoff.
  - Set your FFT grids manually, and choose the one that is used per default for the largest volume (obviously a laborious solution).