diff git a/doc/devel/wavefunctions.rst b/doc/devel/wavefunctions.rst
index 2c5c0ff89ea40dc546435908d54218048f7d78e7..6c1ae0972976995d717656eca21c35bbda992862 100644
 a/doc/devel/wavefunctions.rst
+++ b/doc/devel/wavefunctions.rst
@@ 7,6 +7,12 @@ Wave functions
.. autoclass:: gpaw.wavefunctions.fd.FDWaveFunctions
:members:
+.. autoclass:: gpaw.wavefunctions.pw.PWWaveFunctions
+ :members:
+
+.. autoclass:: gpaw.wavefunctions.pw.PW
+ :members:
+
.. autoclass:: gpaw.wavefunctions.lcao.LCAOWaveFunctions
:members:
@@ 15,4 +21,3 @@ Wave functions
.. autoclass:: gpaw.eigensolvers.eigensolver.Eigensolver
:members:

diff git a/doc/exercises/stress/con_pw.py b/doc/exercises/stress/con_pw.py
index b98fe59d52971e9b209d9f30adaf34b331ee8c55..e320597ba48c994a940e12fb12a2aa707104f96d 100644
 a/doc/exercises/stress/con_pw.py
+++ b/doc/exercises/stress/con_pw.py
@@ 1,12 +1,13 @@
from __future__ import print_function
from ase import Atoms
from ase.build import bulk
from gpaw import GPAW
from gpaw import PW
+from gpaw import GPAW, PW
cell = bulk('Si', 'fcc', a=5.421).get_cell()
a = Atoms('Si2', cell=cell, pbc=True,
 scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))
+a = 5.421
+si = bulk('Si', 'fcc', a=a)
+# or equivalently:
+# b = a / 2
+# from ase import Atoms
+# si = Atoms('Si2', cell=[[0, b, b], [b, 0, b], [b, b, 0]], pbc=True,
+# scaled_positions=[[0, 0, 0], [0.25, 0.25, 0.25]])
for x in [100, 200, 300, 400, 500, 600, 700, 800]:
# for x in [0.24, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.1]:
@@ 16,6 +17,6 @@ for x in [100, 200, 300, 400, 500, 600, 700, 800]:
kpts=(4, 4, 4),
txt='convergence_%s.txt' % x)
 a.set_calculator(calc)
+ si.calc = calc
 print(x, a.get_potential_energy())
+ print(x, si.get_potential_energy())
diff git a/doc/exercises/stress/stress.py b/doc/exercises/stress/stress.py
index cce6b705d3a5791858daee551406c34e41121dbd..cedd6845037e0fab7f027842ba573232fee69760 100644
 a/doc/exercises/stress/stress.py
+++ b/doc/exercises/stress/stress.py
@@ 1,28 +1,22 @@
from __future__ import print_function
from ase import Atoms
+import numpy as np
from ase.build import bulk
from ase.optimize.bfgs import BFGS
from ase.constraints import UnitCellFilter
from gpaw import GPAW
from gpaw import PW
import numpy as np
cell = bulk('Si', 'fcc', a=6.0).get_cell()
+si = bulk('Si', 'fcc', a=6.0)
# Experimental Lattice constant is a=5.421 A
a = Atoms('Si2', cell=cell, pbc=True,
 scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))
calc = GPAW(xc='PBE',
 mode=PW(400),
 # mode=PW(400, cell=a.get_cell()), # fix no of planewaves!
 kpts=(4, 4, 4),
 # convergence={'eigenstates': 1.e10}, # converge tightly!
 txt='stress.txt')
a.set_calculator(calc)
+si.calc = GPAW(xc='PBE',
+ mode=PW(400, dedecut='estimate'),
+ kpts=(4, 4, 4),
+ # convergence={'eigenstates': 1.e10}, # converge tightly!
+ txt='stress.txt')
uf = UnitCellFilter(a)
+uf = UnitCellFilter(si)
relax = BFGS(uf)
relax.run(fmax=0.05) # Consider much tighter fmax!
a = np.dot(a.get_cell()[0], a.get_cell()[0])**0.5 * 2**0.5
print('Relaxed lattice parameter: a = %s A' % a)
+a = np.linalg.norm(si.cell[0]) * 2**0.5
+print('Relaxed lattice parameter: a = {} Ang'.format(a))
diff git a/doc/exercises/stress/stress.rst b/doc/exercises/stress/stress.rst
index 95c386c38163db0ee3896173e08cbc05c1a7a1d3..f964b8b185e8b8d083ea0205a79a8caf3a188ac9 100644
 a/doc/exercises/stress/stress.rst
+++ b/doc/exercises/stress/stress.rst
@@ 1,38 +1,69 @@
.. _stress:
==========================================================
+=================================
Plane wave mode and Stress tensor
==========================================================
+=================================
+
+The major advantage of running DFT calculations on a real space grid, is a
+very efficient parallelization scheme when dealing with large systems.
+However, for small systems it is often faster to use a plane wave basis set
+instead. In this case all quantities are represented by their Fourier
+transforms on the periodic super cell and periodic boundary conditions are
+required. In grid mode, the key convergence parameter is the grid spacing
+`h`, whereas in planewave mode the corresponding parameter is
+`E_{cut}=G_{cut}^2/2` (in atomic units). `G_{cut}` determines the maximum
+size of reciprocal lattice vectors to be included in the plane wave
+expansion.
The major advantage of running DFT calculations on a real space grid, is a very efficient parallelization scheme when dealing with large systems. However, for small systems it is often faster to use a plane wave basis set instead. In this case all quantities are represented by their Fourier transforms on the periodic super cell and periodic boundary conditions are required. In grid mode, the key convergence parameter is the grid spacing `h`, whereas in planewave mode the corresponding parameter is `E_{cut}=G_{cut}^2/2` (in atomic units). `G_{cut}` determines the maximum size of reciprocal lattice vectors to be included in the plane wave expansion.
Converging the plane wave cutoff

+
As a simple example we perform a calculation on bulk Si using the plane wave basis. The following script tests the convergence of the total energy with respect to the cutoff parameter `E_{cut}`.
+As a simple example we perform a calculation on bulk Si using the plane wave
+basis. The following script tests the convergence of the total energy with
+respect to the cutoff parameter `E_{cut}`.
.. literalinclude:: con_pw.py
Note that we use a rather rough kpoint sampling and the default Fermi smearing of 0.1 eV. These parameters should of course be converged, but for now we will just keep them fixed in order to speed up the calculations. At `E_{cut}` = 400 eV, the total energy is seen to be converged to within 5 meV. Now edit the script such that the calculation is performed in grid mode with various values of the grid spacing (comment out the lines :samp:`for x in [100, ...` and :samp:`mode=PW(x)` and uncomment the lines :samp:`for h in [0.24, ...` and :samp:`h=x`).
+Note that we use a rather rough kpoint sampling and the default Fermi
+smearing of 0.1 eV. These parameters should of course be converged, but for
+now we will just keep them fixed in order to speed up the calculations. At
+`E_{cut}` = 400 eV, the total energy is seen to be converged to within 5 meV.
+Now edit the script such that the calculation is performed in grid mode with
+various values of the grid spacing (comment out the lines :samp:`for x in
+[100, ...` and :samp:`mode=PW(x)` and uncomment the lines :samp:`for h in
+[0.24, ...` and :samp:`h=x`).
+
+What grid spacing is needed in order to converge the total energy to within 5
+meV? What is the timing of this calculation compared to the one at 400 eV in
+plane wave mode? Compare the number of grid points and plane waves used in
+these two calculations respectively (look in the output file
+'convergence_400.txt' and the corresponding one for the grid calculation.)
+You can see how long a calculation took by looking at the times recorded in
+the .txt file.
What grid spacing is needed in order to converge the total energy to within 5 meV? What is the timing of this calculation compared to the one at 400 eV in plane wave mode? Compare the number of grid points and plane waves used in these two calculations respectively (look in the output file 'convergence_400.txt' and the corresponding one for the grid calculation.) You can see how long a calculation took by looking at the times recorded in the .txt file.
Optimizing the unit cell

**Warning**: due to difficulties in optimizing cell and positions simultaneously
:class:`ase.constraints.UnitCellFilter` may produce
incorrect results. Always verify obtained structures by means of
performing separate cell
(see :class:`ase.constraints.StrainFilter`)
and positions optimizations (see :mod:`ase.optimize`).
Consider much more tighter fmax than the one used in this tutorial!
+**Warning**: due to difficulties in optimizing cell and positions
+simultaneously :class:`ase.constraints.UnitCellFilter` may produce incorrect
+results. Always verify obtained structures by means of performing separate
+cell (see :class:`ase.constraints.StrainFilter`) and positions optimizations
+(see :mod:`ase.optimize`). Consider a much more tight fmax than the one used
+in this tutorial!
In the :ref:`aluminium_exercise` exercise the lattice constant of bulk Al was found by calculating the total energy at various lattice distances. A nice feature e of the plane wave mode is that it allows a simple implementation of the stress tensor, which can be used to optimize unit unit cells of periodic systems directly. The following script performs such an optimization for bulk Si.
+In the :ref:`aluminium_exercise` exercise the lattice constant of bulk Al was
+found by calculating the total energy at various lattice distances. A nice
+feature e of the plane wave mode is that it allows a simple implementation of
+the stress tensor, which can be used to optimize unit unit cells of periodic
+systems directly. The following script performs such an optimization for bulk
+Si.
.. literalinclude:: stress.py
The calculation uses 12 iterations to find the optimal lattice constant and the relaxation can be viewed with the command line tool :program:`ase`:
+The calculation uses 12 iterations to find the optimal lattice constant and
+the relaxation can be viewed with the command line tool :program:`ase`:
.. highlight:: bash
@@ 40,4 +71,8 @@ The calculation uses 12 iterations to find the optimal lattice constant and the
$ ase gui stress.txt
Since we know the experimental lattice constant, we could probably have calculated the PBE lattice constant faster by fitting a parabola to five points in the vicinity of the expermental lattice constant. However, for complicated unit cells with more than one lattice parameter, the stress tensor becomes a highly valuable tool.
+Since we know the experimental lattice constant, we could probably have
+calculated the PBE lattice constant faster by fitting a parabola to five
+points in the vicinity of the expermental lattice constant. However, for
+complicated unit cells with more than one lattice parameter, the stress
+tensor becomes a highly valuable tool.
diff git a/gpaw/wavefunctions/pw.py b/gpaw/wavefunctions/pw.py
index 98525e98a59489396d5ac998b883d801f172166c..6fe1fef51a9eaa2c8ade01fd8f93527ff5b8a65a 100644
 a/gpaw/wavefunctions/pw.py
+++ b/gpaw/wavefunctions/pw.py
@@ 37,7 +37,7 @@ class PW(Mode):
ecut: float
Planewave cutoff in eV.
 dedecut: float or None
+ dedecut: float or None or 'estimate'
Estimate of derivative of total energy with respect to
planewave cutoff. Used to calculate pulay_stress.
pulay_stress: float or None