""" Make a bandstructure plot using the Harris functional. 
"""
import sys
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
import os

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 Cu 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


bulk = ListOfAtoms([Atom('Cu', (0,    0,     0))] )
a0_fcc = 3.61
b = a0_fcc/2.
bulk.SetUnitCell([(0, b, b), (b, 0, b), (b, b, 0)])

calc = Dacapo(
              kpts=(10,10,10),           # set the k-points Monkhorst-Pack 
              planewavecutoff=340,       # planewavecutoff in eV
              nbands=8,                  # set the number of electronic bands
              usesymm=True ,             # use symmetry to reduce the k-point set
              out='Cu.nc',               # define the out netcdf file
              txtout='Cu.txt' )

bulk.SetCalculator(calc)

calc.StayAliveOff()
calc.SetDensityCutoff(400)

energy = calc.GetPotentialEnergy()
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=6,                  # set the number of electronic bands
              usesymm=True ,             # use symmetry to reduce the k-point set
              out='Cu-Harris.nc',        # define the out netcdf file
              txtout='Cu-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))]

# using the fermilevel from uniform k-point sampling
harrisplot = plotbands(kpoints,eigenvalues,fermilevel)