from ase.build import molecule
from ase.io import write
from gpaw import GPAW

atoms = molecule('H2O')
atoms.pbc = 1
atoms.center(vacuum=3.0)
calc = GPAW(mode='lcao', basis='dzp',
            kpts=[4, 4, 4])

atoms.calc = calc
atoms.get_potential_energy()

nbands = calc.get_number_of_bands()
kpts = calc.get_ibz_k_points()

for k, phase in enumerate(kpts):
    for n in range(nbands):
        psi = calc.get_pseudo_wave_function(band=n, kpt=k)
        # Here you can do /anything/ with psi; it's a numpy array.
        #
        # Since you have complex wavefunctions you may need to
        # do two calls here: data=psi.real and data=psi.imag
        #
        # This only saves the absolute value:
        write('output.kpt{:03d}.band{:03d}.cube'.format(n, k),
              atoms, data=psi)