#!/usr/bin/env python
from Dacapo import Dacapo
from ASE import Atom,ListOfAtoms
from ASE.Visualization.gnuplot import gnuplot
from ASE.Utilities.List import List
import Numeric as num

#Getting ListOfAtoms from old .nc file
bulk=Dacapo.ReadAtoms('Al-fcc.nc')
calc=bulk.GetCalculator()
fermilevel = calc.GetFermilevel()

# setup the Harris calculation

# setup the line from the Gamme point (0,0,0) to the X point (0.5,0.5,0) 
N = 80
Nf = float(N)
kpts = [(x/Nf,x/Nf,0) for x in range(N/2)]
harris = Dacapo(
              kpts=kpts,                 # set the k-points along one direction
              planewavecutoff=340,       # planewavecutoff in eV
              nbands=12,                  # set the number of electronic bands
              usesymm=True ,             # use symmetry to reduce the k-point set
              out='Al-Harris.nc',        # define the out netcdf file
              txtout='Al-Harris.txt',
              infile = 'in.nc' )

bulk.SetCalculator(harris)

harris.StayAliveOff()
harris.SetChargeMixing(False)

# set the density found using the uniform k-point sampling
harris.SetDensityArray(calc.GetDensityArray())

energy1 = harris.GetPotentialEnergy()
kpoints     = harris.GetIBZKPoints()
eigenvalues = [harris.GetEigenvalues(kpt=kpt) for kpt in range(len(kpoints))]

#Plotting the band diagram
def plotbands(kpoints,eigenvalues,fermilevel):
        # make a plot of the eigenvalues
     data = List()
     for kptno in range(len(kpoints)):
         kpointarray = num.array(kpoints[kptno])
         kpoint = num.sqrt(num.sum(kpointarray**2))
         for eigen in eigenvalues[kptno]:
              data.append((kpoint,eigen-fermilevel))
     data.title  = 'Eigenvalues for fcc Al plotted along the Gamma-X point'
     data.xlabel = 'Absolute value of k-point (in units the the reciproval lattice vectors)'
     data.ylabel = 'Eigen values relative to the Fermi energy (eV)'
     data.with   = 'points 3 3'
     bandplot = gnuplot(data)
     return bandplot

# using the fermilevel from uniform k-point sampling

harrisplot = plotbands(kpoints,eigenvalues,fermilevel)