1 #!/usr/bin/env python
   2 from Dacapo import Dacapo
   3 from ASE import Atom,ListOfAtoms
   4 from ASE.Visualization.gnuplot import gnuplot
   5 from ASE.Utilities.List import List
   6 import Numeric as num
   7 
   8 #Getting ListOfAtoms from old .nc file
   9 bulk=Dacapo.ReadAtoms('Al-fcc.nc')
  10 calc=bulk.GetCalculator()
  11 fermilevel = calc.GetFermilevel()
  12 
  13 # setup the Harris calculation
  14 
  15 # setup the line from the Gamme point (0,0,0) to the X point (0.5,0.5,0) 
  16 N = 80
  17 Nf = float(N)
  18 kpts = [(x/Nf,x/Nf,0) for x in range(N/2)]
  19 harris = Dacapo(
  20               kpts=kpts,                 # set the k-points along one direction
  21               planewavecutoff=340,       # planewavecutoff in eV
  22               nbands=12,                  # set the number of electronic bands
  23               usesymm=True ,             # use symmetry to reduce the k-point set
  24               out='Al-Harris.nc',        # define the out netcdf file
  25               txtout='Al-Harris.txt',
  26               infile = 'in.nc' )
  27 
  28 bulk.SetCalculator(harris)
  29 
  30 harris.StayAliveOff()
  31 harris.SetChargeMixing(False)
  32 
  33 # set the density found using the uniform k-point sampling
  34 harris.SetDensityArray(calc.GetDensityArray())
  35 
  36 energy1 = harris.GetPotentialEnergy()
  37 kpoints     = harris.GetIBZKPoints()
  38 eigenvalues = [harris.GetEigenvalues(kpt=kpt) for kpt in range(len(kpoints))]
  39 
  40 #Plotting the band diagram
  41 def plotbands(kpoints,eigenvalues,fermilevel):
  42         # make a plot of the eigenvalues
  43      data = List()
  44      for kptno in range(len(kpoints)):
  45          kpointarray = num.array(kpoints[kptno])
  46          kpoint = num.sqrt(num.sum(kpointarray**2))
  47          for eigen in eigenvalues[kptno]:
  48               data.append((kpoint,eigen-fermilevel))
  49      data.title  = 'Eigenvalues for fcc Al plotted along the Gamma-X point'
  50      data.xlabel = 'Absolute value of k-point (in units the the reciproval lattice vectors)'
  51      data.ylabel = 'Eigen values relative to the Fermi energy (eV)'
  52      data.with   = 'points 3 3'
  53      bandplot = gnuplot(data)
  54      return bandplot
  55 
  56 # using the fermilevel from uniform k-point sampling
  57 
  58 harrisplot = plotbands(kpoints,eigenvalues,fermilevel)