If you are doing energy-volume calculations or cell shape and volume relaxations you must understand the Pulay stress, and related problems.
The Pulay stress arises from the fact that the plane wave basis set is not complete with respect to changes of the volume. Thus, unless absolute convergence with respect to the basis set has been achieved - the diagonal components of the stress tensor are incorrect. This error is often called ``Pulay stress''. The error is almost isotropic (i.e. the same for each diagonal component), and for a finite basis set it tends to decrease volume compared to fully converged calculations (or calculations with a constant energy cutoff).
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 cutoff region becomes an elipsoide. Restarting VASP after a volume relaxation causes VASP to adopt a new ``spherical'' cutoff sphere and thus the energy changes discontinuously (see section 6.14).
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, see section 6.14), 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. But calculations with many plane wave codes have shown that such calculations give very reliable results for the lattice constant and the bulk modulus and other elastic properties even at relatively small energy cutoffs. Constant energy cut-off calculations are less prone to errors cause 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 in the limit of very small distortions the energy changes obtained with this method are equivalent with that 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.