1 """ workfunction and dipole correction for a Al slab 
   2 """
   3 
   4 import os
   5 import Numeric as num
   6 from ASE.Visualization.gnuplot import gnuplot
   7 from ASE.Visualization.VTK import VTKPlotArray
   8 from Dacapo import Dacapo
   9 from ASE import Atom, ListOfAtoms
  10 from ASE.Utilities.List import List
  11 import os.path,sys
  12 from ASE.IO.Cube import WriteCube
  13 
  14 sys.path.remove('/home/camp/lhansen/CamposASE')
  15 sys.path.remove('/usr/lib/python2.2/site-packages/CamposASE')
  16 
  17 # from Visualization import *
  18 
  19 # check if we have already calculated the nc file
  20 if not os.path.isfile('Al111-O.nc'):
  21 
  22     # setup the Al slab
  23     alslab =  ListOfAtoms([ Atom('Al',(0, 0, 0)),
  24                         Atom('Al',(0, 0.5, 0)),
  25                         Atom('Al',(0.5, 0, 0)),
  26                         Atom('Al',(0.5, 0.5, 0)),
  27                         Atom('Al',(0.3333333, 0.166666667, 0.2338/1.25)),
  28                         Atom('Al',(0.3333333, 0.666666667, 0.2338/1.25)),
  29                         Atom('Al',(0.8333333, 0.166666667, 0.2338/1.25)),
  30                         Atom('Al',(0.8333333, 0.666666667, 0.2338/1.25 ))])
  31 
  32     alslab.append(Atom('O',(0.0,0.0, -0.12),tag=1))
  33 
  34     unitcell = [[5.72756492761103, 0, 0],
  35                 [-2.86378246380552, 4.96021672913593, 0],
  36                 [0,0,12.5]]
  37 
  38     alslab.SetUnitCell(unitcell)
  39 
  40     calc = Dacapo(planewavecutoff = 340,
  41                   densitycutoff = 500,
  42                   dipole = True,
  43                   nbands = 30, 
  44                   out = 'Al111-O.nc', 
  45                   txtout = 'Al111-O.txt') 
  46 
  47 
  48     alslab.SetCalculator(calc)
  49 
  50 else: 
  51     alslab = Dacapo.ReadAtoms('Al111-O.nc') 
  52     calc = alslab.GetCalculator()
  53 
  54 
  55 electro = calc.GetElectrostaticPotential()
  56 
  57 nx,ny,nz = num.shape(electro) 
  58 
  59 # get the x-axis 
  60 unitcell = alslab.GetUnitCell()
  61 xaxis = [float(z)/float(nz)*unitcell[2,2] for z in range(nz)]
  62 
  63 # xy average of potential
  64 xyaverage = [num.sum(num.sum(electro[:,:,z]))/(nx*ny) for z in range(nz)]
  65 potential = List(zip(xaxis,xyaverage))
  66 potential.legend = 'Electrostatic potential'
  67 potential.ylabel = 'Electrostatic potential (eV)'
  68 potential.xlabel = 'Distance along z-axis (Angstrom)'
  69 
  70 # add fermilevel to the plot
  71 fermilevel = calc.GetFermilevel()
  72 fermilevel = List(zip(xaxis,[fermilevel for z in range(nz)]))
  73 fermilevel.legend = 'Fermi level'
  74 
  75 # make the plots
  76 plot = gnuplot(potential) 
  77 fermiplot = gnuplot(fermilevel,parent=plot) 
  78 plot.Update()
  79 
  80 electro = calc.GetDensityArray()
  81 
  82 # vtkplot = VTKPlotArray(electro,unitcell)
  83 WriteCube(alslab,electro[:,:,:],'test.cube')