Pulay stress: Difference between revisions

From VASP Wiki
VisualWikitext
No edit summary
No edit summary
Line 8: Line 8:
The Pulay stress may be corrected in multiple ways. Generally, by calculating the relaxed structure with a larger basis-set by increasing the {{TAG|ENCUT}}:
The Pulay stress may be corrected in multiple ways. Generally, by calculating the relaxed structure with a larger basis-set by increasing the {{TAG|ENCUT}}:
#Set {{TAG|ENCUT}} to <math>1.3\times</math> the default cutoff, or {{TAG|PREC}}=''High'' in VASP.4.4.
#Set {{TAG|ENCUT}} to <math>1.3\times</math> the default cutoff, or {{TAG|PREC}}=''High'' in VASP.4.4.
#Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.  
#Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.<ref>The most reliable relaxed structures 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.</ref>


One thing which is important to understand is 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 almost uniform and changes the diagonal elements of the stress tensor only by a constant amount.  
One thing which is important to understand is 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 almost uniform and changes the diagonal elements of the stress tensor only by a constant amount.  
Line 28: Line 28:


=== Step 4a. ===
=== Step 4a. ===
If there are still problems, try further increasing the {{TAG|ENCUT}}. N<ref>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.</ref>
If there are still problems, try further increasing the {{TAG|ENCUT}}.  


=== Step 4b. ===
=== Step 4b. ===

Revision as of 14:47, 1 July 2024

Pseudopotential calculations are 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, i.e.g 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 creates stress that decreases the volume compared to fully converged calculations, due to the diagonal components of the stress tensor being incorrect. This is the "Pulay stress".[2] Figures 1 and 2 shows these discontinuities for diamond for prematurely truncated energy cutoffs and k-point meshes, respectively, in comparison to converged curves. It may be corrected for using a few different methods in VASP.

Fig 1. Total energy vs lattice parameter for converged and unconverged plane wave energy cutoffs. Diamond in a primitive cell - 2x2x2 k-point mesh.
Fig 2. Total energy vs lattice parameter for converged and unconverged k-point meshes. Diamond in a primitive cell - 250 eV energy cutoff.

How to correct in VASP

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

  1. Set ENCUT to the default cutoff, or PREC=High in VASP.4.4.
  2. Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.[3]

One thing which is important to understand is 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 almost uniform and changes the diagonal elements of the stress tensor only by a constant amount.

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

Accurate Bulk Relaxations - Structure Relaxation

One way is to very accurately relax the structure in a series of steps.

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, a few more additional steps may be tried:

Step 4a.

If there are still problems, try further increasing the ENCUT.

Step 4b.

If the default energy cutoff must be used for the relaxations, then an alternative procedure may be followed. [5]

  1. Perform a single point calculation for starting structure with a higher ENCUT. In the OUTCAR file, good estimation for the Pulay stress is given as a negative pressure, e.g.:
  external pressure =    -100.29567 kB
  1. Add PSTRESS to the INCAR.
  2. Set it equal to the Pulay stress given by the 'external pressure' in the OUTCAR.

Accurate Bulk Relaxations - Equation of State Fitting

An alternative way to avoid relaxing the volume is to relax the structure for a fixed set of volumes, i.e. multiple POSCAR files. These are then fitted to an equation of state[6], e.g. Murnaghan.[7]. This has the advantage of also providing the bulk modulus, and we have found it may be used safely with the default energy cut-off.

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

  1. Select one cell and relax from starting structure, keeping the volume fixed (ISIF=4; ISMEAR=0 or 1).
  2. Start a second relaxation (structure or electronic 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).
  3. 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).

does this procedure actually make sense?

Possible issues and advice on how to address it

It is very common among new users of plane wave codes to obtain jagged energy vs. volume plots. This is commonly due to two reasons:

  1. 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.[9]
  2. Discontinuity of FFT grids

    Different precisions of the FFT grid defined by PREC may be used, e.g. Normal or Accurate.[10] 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.

    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

  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. ↑ The Pulay stress and related problems affect the behavior of VASP and any plane wave code in several besides the stress tensor. All volume/cell shape relaxation algorithms implemented in VASP work with a constant basis set. In that way, all energy changes are strictly consistent with the calculated stress tensor, and this, in turn, results in an underestimation of the equilibrium volume unless a large plane wave cutoff is used. Keeping the basis set constant during relaxations has also a strange effect on the basis set. Initially, all G-vectors within a sphere are included in the basis. If the cell shape relaxation is performed, the direct and reciprocal lattice vectors change. This means that although the number of reciprocal G-vectors in the basis is kept fixed, the length of the G-vectors changes, changing indirectly the energy cutoff. Or to be more precise, the shape of the cutoff region changes from a sphere to an ellipsoid. Restarting VASP after a volume relaxation causes VASP to adopt a new "spherical" cutoff sphere and thus the energy changes discontinuously.
  3. ↑ The most reliable relaxed structures 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.
  4. ↑ Calculations with many plane-wave codes have shown that such calculations yield reliable results for the lattice constants and the bulk modulus and other elastic properties even at relatively modest energy cutoffs. Constant energy cut-off calculations are less prone to errors caused by the basis set incompleteness than constant basis set calculations. But it should be kept in mind that volume and cell shape changes must be rather large in order to obtain reliable results from this method because within the limit of very small distortions, the energy changes obtained with this method are equivalent to those obtained from the stress tensor and are therefore affected by the Pulay stress. Only volume changes of the order of 5-10 % guarantee that the errors introduced by the basis set incompleteness are averaged out.
  5. ↑ 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.
  6. ↑ 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.
  7. ↑ Murnaghan Equation of State, www.wikipedia.org (2024)
  8. ↑ 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.
  9. ↑ This solves the problem, because the criterion for including a plane wave in the basis set is . 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.
  10. ↑ For PREC=Accurate, the FFT grids are chosen such that 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 .
Retrieved from "https://www.vasp.at/wiki/index.php?title=Pulay_stress&oldid=25488"