from ase.data.molecules import *
m = molecule('BeH')
from gpaw import GPAW
calc=GPAW(idiotproof=False)
m.center(2.5)
m.set_calculator(calc)
calc.initialize(m)

import numpy as np
f_si = calc.wfs.setups[0].calculate_initial_occupation_numbers(magmom=0,
                                                               hund=0,
                                                               charge=0,
                                                               nspins=2
                                                               )
print f_si
def occupations(*args, **kwargs):
    return f_si

calc.wfs.setups[0].calculate_initial_occupation_numbers = occupations
m.get_potential_energy()