1 """ Make a bandstructure plot using the Harris functional. 
   2 """
   3 import sys
   4 from Dacapo import Dacapo
   5 from ASE import Atom,ListOfAtoms
   6 from ASE.Visualization.gnuplot import gnuplot
   7 from ASE.Utilities.List import List
   8 import Numeric as num
   9 import os
  10 
  11 def plotbands(kpoints,eigenvalues,fermilevel): 
  12 	# make a plot of the eigenvalues
  13 
  14      data = List()
  15      for kptno in range(len(kpoints)):
  16 
  17          kpointarray = num.array(kpoints[kptno])
  18 	 kpoint = num.sqrt(num.sum(kpointarray**2))
  19 	    
  20 	 for eigen in eigenvalues[kptno]:
  21 	      data.append((kpoint,eigen-fermilevel))
  22 
  23 
  24      data.title  = 'Eigenvalues for fcc Cu plotted along the Gamma-X point'
  25      data.xlabel = 'Absolute value of k-point (in units the the reciproval lattice vectors)'
  26      data.ylabel = 'Eigen values relative to the Fermi energy (eV)'
  27      data.with   = 'points 3 3'
  28      bandplot = gnuplot(data)
  29      return bandplot
  30 
  31 
  32 bulk = ListOfAtoms([Atom('Cu', (0,    0,     0))] )
  33 a0_fcc = 3.61
  34 b = a0_fcc/2.
  35 bulk.SetUnitCell([(0, b, b), (b, 0, b), (b, b, 0)])
  36 
  37 calc = Dacapo(
  38               kpts=(10,10,10),           # set the k-points Monkhorst-Pack 
  39               planewavecutoff=340,       # planewavecutoff in eV
  40               nbands=8,                  # set the number of electronic bands
  41               usesymm=True ,             # use symmetry to reduce the k-point set
  42               out='Cu.nc',               # define the out netcdf file
  43               txtout='Cu.txt' )
  44 
  45 bulk.SetCalculator(calc)
  46 
  47 calc.StayAliveOff()
  48 calc.SetDensityCutoff(400)
  49 
  50 energy = calc.GetPotentialEnergy()
  51 fermilevel = calc.GetFermilevel()
  52 
  53 # setup the Harris calculation
  54 
  55 # setup the line from the Gamme point (0,0,0) to the X point (0.5,0.5,0) 
  56 N = 80
  57 Nf = float(N) 
  58 kpts = [(x/Nf,x/Nf,0) for x in range(N/2)]
  59 harris = Dacapo(
  60               kpts=kpts,                 # set the k-points along one direction
  61               planewavecutoff=340,       # planewavecutoff in eV
  62               nbands=6,                  # set the number of electronic bands
  63               usesymm=True ,             # use symmetry to reduce the k-point set
  64               out='Cu-Harris.nc',        # define the out netcdf file
  65               txtout='Cu-Harris.txt', 
  66               infile = 'in.nc' ) 
  67 
  68 bulk.SetCalculator(harris) 
  69 
  70 harris.StayAliveOff()
  71 harris.SetChargeMixing(False) 
  72 
  73 # set the density found using the uniform k-point sampling
  74 harris.SetDensityArray(calc.GetDensityArray())
  75 
  76 energy1 = harris.GetPotentialEnergy()
  77 kpoints     = harris.GetIBZKPoints()
  78 eigenvalues = [harris.GetEigenvalues(kpt=kpt) for kpt in range(len(kpoints))]
  79 
  80 # using the fermilevel from uniform k-point sampling
  81 harrisplot = plotbands(kpoints,eigenvalues,fermilevel)