Latest revision as of 09:11, 30 August 2024

Pulay stress is unphysical stress resulting from unconverged calculations with respect to the basis set. It distorts the cell structure, decreasing it from the equilibrium volume and creating difficulties in volume relaxation. The resultant energy vs. volume curves, cf. Figure 1 (top), are jagged and special care must be taken to obtain reasonable structures, cf. Volume relaxation. In this article, the computational origin of this is discussed.

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.

Introduction

The energy for a periodic system, e.g. band structures, is calculated using a finite number of plane waves and a finite number of k-points. A fixed number of plane waves or plane wave energy cutoff may be used to set a constant basis.^[1] In VASP, a constant energy cutoff is used, cf. ENCUT. The number of plane waves N_PW (Note: the number of plane waves in VASP can be found using by searching for NPLWV in the OUTCAR file) is related to the energy cutoff E_cutoff and the size of the cell Ω₀:

N_{PW} \propto\ \Omega_0\ E_{cutoff}^{3/2}

N_PW is constant in a relaxation calculation, which means that E_cutoff must change to compensate for changes in Ω₀. All the initial G-vectors within a sphere are included in the basis. However, when comparing cells of different sizes, i.e. during a relaxation, the cell shape is relaxed, so 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. This can 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.^[2] All energy changes are strictly consistent with the stress tensor; 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 equilibrium volume (cf. Fig. 1 (top)), due to the diagonal components of the stress tensor being incorrect. This is called the Pulay stress.

The pressure of the cell, being proportional to the trace of the stress tensor, can be used to visualize this. When the cell volume is below the equilibrium volume, the pressure is positive; contrastingly, it is negative when above the equilibrium volume, so at equilibrium, this is zero. Plotting the magnitude of the pressure vs. volume curve and the total energy allows comparison between these two minima. In Figure 1 it is clear that the the absolute pressure-volume and energy-volume minima coincide for a converged basis, while the pressure equilibrium is much lower than the energy equilibrium for the unconverged basis. This is the effect of the Pulay stress.

Further explanation

As mentioned previously, N_PW is constant in a relaxation calculation, which means that E_cutoff must change to compensate for changes in Ω₀. This is illustrated in Fig. 2. The initial G-vectors within a sphere are included within the basis.

When the cell volume increases (V₁ < V₁), the number of G-vectors in reciprocal space remains constant, but their length increases (cf. Fig. 2 (top)). This effectively results in a change of basis, leading to (E_{cutoff, 1} > E_{cutoff, 2}). This basis remains constant for the duration of the relaxation. However, if the calculation is then restarted, the basis is reset. This means that the number of G-vectors is greater for the larger, real-space cell. One effect of this is that there are more real-space grid points. However, the corresponding reciprocal space decreases.

Contrastingly, see Fig. 2 (bottom), when the volume decreases on relaxation (V₁ > V₁), the length of the G-vectors decreases. The effective E_cutoff should increase but this does not improve the situation, as it creates an artificial pressure. The reciprocal space grid points are effectively sparser. If the calculation restarts, the basis is reset, so the number of G-vectors decreases for the smaller real space cell.

Alternatively, the shape of the cell could change. As the shape changes, the G-vectors continue to be directed along the lattice coordinates, meaning that some shorten while others lengthen, see Fig. 3. This results in a shift from a spherical basis, where all G-vectors are of equal length, to one where some are stretched and others compressed, i.e. an ellipsoid. This changes the effective E_cutoff along each lattice parameter. On resetting the calculation, the cutoff is once again spherical. This draws an analogy to the symmetry breaking of the Bravais lattice seen for gradient-corrected functionals (cf. GGA_COMPAT), where the spherical symmetry of the G-vectors is broken for non-cubic cells.

References

[gomesdacosta:nielsen:kunc:1986-1] P. Gomes Dacosta, O. Nielsen, and K. Kunc, Stress theorem in the determination of static equilibrium by the density functional method, J. Phys. C: Solid State Phys. 19, 3163 (1986).

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

[1]

[2]

@@ Line 1: / Line 1: @@
-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.{{cite|payne:francis:1990}} 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".<ref>The Pulay stress and related problems affect the behavior of VASP and any plane wave code in several ways: First, it evidently affects the stress tensor calculated by VASP, i.e. the diagonal components of the stress tensor are incorrect, unless the energy cutoff is very large (ENMAX=1.3 *default is usually a safe setting to obtain a reliable stress tensor). In addition, it should be noted that 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 some strange effect on the basis set. Initially, all G-vectors within a sphere are included in the basis. If the cell shape relaxation starts 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 becomes an ellipsoid. Restarting VASP after a volume relaxation causes VASP to adopt a new "spherical" cutoff sphere and thus the energy changes discontinuously.
+Pulay stress is unphysical stress resulting from unconverged calculations with respect to the basis set. It distorts the cell structure, decreasing it from the equilibrium volume and creating difficulties in volume relaxation. The resultant energy vs. volume curves, cf. Figure 1 (top), are jagged and special care must be taken to obtain reasonable structures, cf. [[volume relaxation|Volume relaxation]]. In this article, the computational origin of this is discussed. [[File:Pressure_energy_volume.png|400px|thumb|Figure 1. Total energy (left y-axis) and absolute pressure (right y-axis) vs. lattice parameter. Equilibrium lattice parameters for energy and pressure are shown. These coincide when Pulay stress is eliminated. ENCUT = 250 eV (top - unconverged) and 540 eV (bottom - converged). Diamond in a primitive cell - 2x2x2 k-point mesh.]]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.
-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 it, therefore, changes the diagonal elements of the stress tensor only by a certain constant amount (see below). In addition, many calculations have shown that Pulay stress-related problems can also be reduced by performing calculations at different volumes using the same energy cutoff for each calculation (this is what VASP does by default) and fitting the final energies to an equation of state. This of course implies that the number of basis vectors is different at each volume. 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 changes 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.</ref> 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.
+= Introduction =
+The energy for a periodic system, e.g. band structures, is calculated using a finite number of plane waves and a finite number of k-points. A fixed number of plane waves or plane wave energy cutoff may be used to set a constant basis.{{cite|gomesdacosta:nielsen:kunc:1986}} In VASP, a constant energy cutoff is used, cf. {{TAG|ENCUT}}. The number of plane waves <b><i><span>N</span></b><sub>PW</sub></i> (Note: the number of plane waves in VASP can be found using by searching for ''NPLWV'' in the {{TAG|OUTCAR}} file) is related to the energy cutoff <b><i><span>E</span></b><sub>cutoff</sub></i> and the size of the cell <span>'''&Omega;'''</span><sub>0</sub>:
+::<math> N_{PW} \propto\ \Omega_0\ E_{cutoff}^{3/2} </math>
-[[File:ENCUT_comp.png|300px|thumb|Fig 1. Total energy vs lattice parameter for converged and unconverged plane wave energy cutoffs. Diamond in a primitive cell. 2x2x2 k-point mesh.]]
+<b><i><span>N</span></b><sub>PW</sub></i> is constant in a relaxation calculation, which means that <b><i><span>E</span></b><sub>cutoff</sub></i> must change to compensate for changes in <span>'''&Omega;'''</span><sub>0</sub>. All the initial G-vectors within a sphere are included in the basis. However, when comparing cells of different sizes, i.e. during a relaxation, the cell shape is relaxed, so 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. This can 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.{{cite|payne:francis:1990}} All energy changes are strictly consistent with the stress tensor; 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 equilibrium volume (cf. Fig. 1 (top)), due to the diagonal components of the stress tensor being incorrect. This is called the ''Pulay stress''.
-[[File:Kpoint_comp.png|300px|thumb|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 pressure of the cell, being proportional to the trace of the stress tensor, can be used to visualize this. When the cell volume is below the equilibrium volume, the pressure is positive; contrastingly, it is negative when above the equilibrium volume, so at equilibrium, this is zero. Plotting the magnitude of the pressure vs. volume curve and the total energy allows comparison between these two minima. In Figure 1 it is clear that the the absolute pressure-volume and energy-volume minima coincide for a converged basis, while the pressure equilibrium is much lower than the energy equilibrium for the unconverged basis. This is the effect of the Pulay stress.
-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.
-#Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.
-=== Accurate bulk relaxations with internal parameters (one) ===
+= Further explanation =
+As mentioned previously, <b><i><span>N</span></b><sub>PW</sub></i> is constant in a relaxation calculation, which means that <b><i><span>E</span></b><sub>cutoff</sub></i> must change to compensate for changes in <span>'''&Omega;'''</span><sub>0</sub>. This is illustrated in Fig. 2. The initial G-vectors within a sphere are included within the basis.
-If volume relaxations are necessary, follow this procedure:
+When the cell volume increases (<b><span>V</span></b><sub>1</sub> < <b><span>V</span></b><sub>1</sub>), the number of G-vectors in reciprocal space remains constant, but their length increases (cf. Fig. 2 (top)). This effectively results in a change of basis, leading to (<b><i><span>E</span></b><sub>cutoff, 1</sub></i> > <b><i><span>E</span></b><sub>cutoff, 2</sub></i>). This basis remains constant for the duration of the relaxation. However, if the calculation is then restarted, the basis is reset. This means that the number of G-vectors is greater for the larger, real-space cell. One effect of this is that there are more real-space grid points. However, the corresponding reciprocal space decreases.
-=== Step 1. ===
+Contrastingly, see Fig. 2 (bottom), when the volume decreases on relaxation (<b><span>V</span></b><sub>1</sub> > <b><span>V</span></b><sub>1</sub>), the length of the G-vectors decreases. The effective <b><i><span>E</span></b><sub>cutoff</sub></i> should increase but this does not improve the situation, as it creates an artificial pressure. The reciprocal space grid points are effectively sparser. If the calculation restarts, the basis is reset, so the number of G-vectors decreases for the smaller real space cell.
-Relax from the starting structure with Gaussian or Methfessel-Paxton smearing ({{TAG|ISMEAR}} = 0 or 1).
-=== Step 2. ===
-Copy the {{TAG|CONTCAR}} to {{TAG|POSCAR}} and relax the structure again.
-=== Step 3. ===
-Change the smearing method to the tetrahedron method (i.e. {{TAG|ISMEAR}}=-5) and perform a single point calculation, i.e. no relaxation of structure. These are will give highly accurate energies.
-This should give good energies. If there are still problems, a few more additional steps can be tried.
+[[File:Pulay_stress_grids.png|700px|thumb|centre|Figure 2. Cell shape and lattice positions are kept constant, while the volume <b><span>V</span></b> is free to change (ISIF = 7). The initial volume <b><span>V</span></b><sub>1</sub> changes to the final volume <b><span>V</span></b><sub>2</sub>. Two cases are given, one for volume increasing on relaxation (top) and one for it decreasing (bottom). The change in real space is given on the left, while the change in reciprocal space and the subsequent effect on the G-vectors is given on the right. Blue is the initial basis, while red is the new, restarted basis. The relation between <b><i><span>E</span></b><sub>cutoff</sub></i>, <b><i><span>N</span></b><sub>PW</sub></i>, and G-vectors is given for the initial and final volumes.]]
-=== Step 4a. ===
+Alternatively, the shape of the cell could change. As the shape changes, the G-vectors continue to be directed along the lattice coordinates, meaning that some shorten while others lengthen, see Fig. 3. This results in a shift from a spherical basis, where all G-vectors are of equal length, to one where some are stretched and others compressed, i.e. an ellipsoid. This changes the effective <b><i><span>E</span></b><sub>cutoff</sub></i> along each lattice parameter. On resetting the calculation, the cutoff is once again spherical. This draws an analogy to the symmetry breaking of the Bravais lattice seen for gradient-corrected functionals (cf. {{TAG|GGA_COMPAT}}), where the spherical symmetry of the G-vectors is broken for non-cubic cells.
-If there are still problems, try further increasing the {{TAG|ENCUT}}.<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>
-=== Step 4b. ===
+[[File:Shape_pulay_grids.png|700px|thumb|centre|Figure 3. Cell volume and lattice positions are kept constant, while the shape is free to change (ISIF = 5). The shape changes from cubic to hexagonal. The blue spherical basis changes to the red ellipsoid basis, along the direction of the sheer. On restarting, a spherical basis returns.]]
-If the default energy cutoff must be used for the relaxations, then an alternative procedure may be followed.
-<ref>From now on all ''STRESS'' output of VASP is corrected by simply subtracting {{TAG|PSTRESS}}. In addition, all volume relaxations will take {{TAG|PSTRESS}} into account.
-It should be said again, use {{TAG|PSTRESS}} only if increasing the cutoff is not a viable option for some reason.</ref>
-<ol>
-<li>Perform a single point calculation for starting structure with a higher {{TAG|ENCUT}}. In the {{TAG|OUTCAR}} file, good estimation for the Pulay stress is given as a negative pressure, e.g.:<li>
-   external pressure =    -100.29567 kB
-<li>Add {{TAG|PSTRESS}} to the INCAR.<li>
-<li>Set it equal to the Pulay stress given by the 'external pressure' in the OUTCAR.<li>
-</ol>
-=== Accurate bulk relaxations with internal parameters (two) ===
+[[Category:Ionic minimization]][[Category:Howto]]
-It is possible to avoid volume relaxation in many cases:
-*Relax the structure (cell shape and internal parameters) for a set of fixed volumes ({{TAG|ISIF}}=4).
-*Fit to an equation of state <ref>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.</ref>, e.g. Murnaghan.{{cite|murnaghan:web}} <ref>This has the advantage of also providing the bulk modulus, and we have found in the past that it may be used safely with the defauly cut-off.</ref>
-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 ({{TAG|ISIF}}=4; {{TAG|ISMEAR}}=0 or 1).
-*Start a second relaxation (structure or electronic relaxation?) from the previous {{TAG|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 ({{TAG|ISMEAR}}=-5), to get very accurate unambiguous energies (no relaxation for the final run).
-== 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 <math>\vert {\bf G} + {\bf k} \vert < {\bf G}_{\rm cut}</math>. 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 {{TAG|PREC}}=Normal. For {{TAG|PREC}}=Accurate the FFT grids are chosen such that <math> {\bf H} \vert \phi> </math> is exactly evaluated. For {{TAG|PREC}}=Normal the FFT grids are set to 3/4 of the value that is required for an exact evaluation of <math> {\bf H} \vert \phi> </math>. 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).
-----
-[[Category:Ionic minimization]][[Category:Howto]]
 ==References==