Pulay stress: Difference between revisions

From VASP Wiki
 
(120 intermediate revisions by the same user not shown)
Line 1: Line 1:
Pulay stress is unphysical stress coming from unconverged calculations with respect to basis set. It distorts the cell structure, decreasing it from the equilibrium volume and resulting in jagged energy vs. volume curves, cf. Figures 1 and 2. In this article, the computational origin of this is briefly discussed, and then solutions for correcting the Pulay stress. 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.
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.
 
[[File:ENCUT_comp.png|400px|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.]]
[[File:Kpoint_comp.png|400px|thumb|Fig 2. Total energy vs. lattice parameter for converged and unconverged k-point meshes. Diamond in a primitive cell - 180 eV energy cutoff.]]


= Introduction =  
= 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 is related to the energy cutoff and the size of the cell:
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} ∝  </math>
::<math> N_{PW} \propto\ \Omega_0\ E_{cutoff}^{3/2} </math>
 
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, due to the diagonal components of the stress tensor being incorrect. This is called 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.
 
= 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 {{TAG|ENCUT}} until convergence is reached:
#Set {{TAG|ENCUT}} to <math>1.3\times</math> the default cutoff, or {{TAG|PREC}}=''High'' in VASP.4.4. N.B. this is the lower recommended limit.
#Re-run VASP with the default cutoff to obtain the final relaxed positions and cell parameters.
#Further increase the {{TAG|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 ({{TAG|ISIF}} = 7; {{TAG|IBRION}} = 1, 2, or 3) with Gaussian or Methfessel-Paxton smearing ({{TAG|ISMEAR}} = 0 or 1). N.B. other settings can be used for {{TAG|ISIF}}, e.g.  3, 6, or 8 but it is not recommended to relax more than one degree of freedom at a time.
 
=== 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.
 
=== (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. ===
<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''.  
Try further increasing the {{TAG|ENCUT}}. Alternatively, improve the FFT grid by setting {{TAG|PREC}}=Accurate.


=== Step 4b. ===
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.
To avoid additional computational cost due to increased cutoff energy, the ''STRESS'' output in VASP may be corrected using {{TAG|PSTRESS}}. The Pulay stress is only weakly dependent on volume and ionic configuration; it is mainly determined by the composition and basis set. A good estimation for it is given in the output, e.g.:
  external pressure =    -100.29567 kB
The difference in this pressure (between the desired and a very large {{TAG|ENCUT}}) may then be used to correct for the Pulay stress. {{TAG|PSTRESS}} is set to this difference, then 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 two single point calculations, one for the default and one for the higher {{TAG|ENCUT}}.</li>
<li>Find the external pressure in the {{TAG|OUTCAR}} file, e.g.:</li>
</ol>
  external pressure =    -1311.32 kB
 
  external pressure =    -95.66 kB
<ol type = 'i' start = "3">
<li>Find the difference between these external pressures. This is a good approximation of the Pulay stress. E.g.</li>
</ol>
  difference in pressure = -1215.66 kB
<ol type = 'i' start = "4">
<li>Set {{TAG|PSTRESS}} equal to this difference in the OUTCAR file, i.e.:
</li>
</ol>
  {{TAG|PSTRESS}} = -1215.66
This results in structures similar to the higher cutoff at the cost of the default cutoff. We reiterate that {{TAG|PSTRESS}} should only be used if the higher cutoff is not a viable option, as this only improves the structure and not the energy.


== 2. Equation of state fitting ==
= Further explanation =
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 {{TAG|POSCAR}} files. These are then fitted to an equation of state, e.g. Murnaghan.{{cite|murnaghan:web}} 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.
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.  


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:
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. ===
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 ({{TAG|ISIF}}=4; {{TAG|ISMEAR}}=0 or 1).


=== Step 2. ===
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.  
Copy the {{TAG|CONTCAR}} to {{TAG|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. {{TAG|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 =
[[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.]]


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:
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.


<ol>
[[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.]]
<li>'''Basis set incompleteness'''</p><p>
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.
<ol type='a'>
<li>Use a larger plane wave cutoff, cf. Fig. 1. This is usually the preferred and simplest solution.</li>
<li>Use more k-points, cf. Fig. 2.<ref>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.</ref></li></ol>
</li>


<li>'''Discontinuity of FFT grids'''</p><p>
[[Category:Ionic minimization]][[Category:Howto]]
Different precisions of the FFT grid defined by {{TAG|PREC}} may be used, e.g. Normal or Accurate.<ref>For {{TAG|PREC}}=Accurate, the FFT grids are chosen such that <math> {\bf H} \vert \phi> </math> is exactly evaluated. Whereas, 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>.</ref> 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 {{cite|vasp:intro:lecture:web}}.
<ol type='a'>
<li>Use {{TAG|PREC}}=Accurate</li>
<li>Increase the plane wave cutoff.</li>
<li>Set your FFT grids manually, and choose the one that is used per default for the largest volume.</li></ol>
</li>
</ol>
<!--[[Category:Ionic minimization]][[Category:Howto]]-->


==References==
==References==

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.

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.

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 NPW (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 Ecutoff and the size of the cell Ω0:

NPW is constant in a relaxation calculation, which means that Ecutoff must change to compensate for changes in Ω0. 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, NPW is constant in a relaxation calculation, which means that Ecutoff must change to compensate for changes in Ω0. This is illustrated in Fig. 2. The initial G-vectors within a sphere are included within the basis.

When the cell volume increases (V1 < V1), 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 (Ecutoff, 1 > Ecutoff, 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 (V1 > V1), the length of the G-vectors decreases. The effective Ecutoff 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.

Figure 2. Cell shape and lattice positions are kept constant, while the volume V is free to change (ISIF = 7). The initial volume V1 changes to the final volume V2. 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 Ecutoff, NPW, and G-vectors is given for the initial and final volumes.

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 Ecutoff 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.

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.

References