import Gnuplot as gp
import Numeric as num
from Dacapo import Dacapo
from BuildFCC import FCC111
from ASE import Atom


a = 4.05
h = 12.0
fcc = FCC111('Al', a, 2, h)

# Add the adsorbate:
fcc.append(Atom('H', (0, 0, 1.55)))

calc = Dacapo(nbands=2 * 5,
              kpts=(4, 4, 1),
              dipole=True,
              planewavecutoff=340,
              densitycutoff=400)
fcc.SetCalculator(calc)

v = calc.GetElectrostaticPotential()
nx, ny, nz = v.shape
vz = num.sum(num.sum(v)) / (nx * ny)
z = num.arange(nz) * h / nz
plot = gp.Gnuplot(persist=True)
plot.plot(gp.Data(z, vz, with='lines lw 3'),
          gp.Func(str(calc.GetFermilevel()), with='lines lw 3'))