from gpaw import GPAW
from ase.data.molecules import molecule
from gpaw.utilities.tools import tri2full

calc = GPAW(mode='lcao', basis='dzp')
system = molecule('H2O')
system.center(vacuum=2.5)

system.set_calculator(calc)
system.get_potential_energy()

H_MM = calc.wfs.eigensolver.calculate_hamiltonian_matrix(calc.hamiltonian,
                                                         calc.wfs,
                                                         calc.wfs.kpt_u[0])
tri2full(H_MM)

index_of_oxygen_atom = 0 # or something

M = calc.wfs.basis_functions.M_a[index_of_oxygen_atom]
nao = calc.wfs.setups[index_of_oxygen_atom].niAO

interesting_hamiltonian_slice = H_MM[M:M + nao, :]

print interesting_hamiltonian_slice