from gpaw import GPAW
from ase.data.molecules import molecule

if 0:
    system = molecule('H2O')
    system.center(vacuum=3.0)

    calc = GPAW()

    system.set_calculator(calc)
    system.get_potential_energy()

    calc.write('h2o.gpw')

calc = GPAW('h2o.gpw')

vxc_g = calc.hamiltonian.finegd.zeros()
calc.density.interpolate()
nt_g = calc.density.nt_sg[0]

Exc = calc.hamiltonian.xc.get_energy_and_potential(nt_g, vxc_g)

print Exc
print nt_g[7, 7, 4:7]