1 #!/usr/bin/env python
   2 # =======================================
   3 #    STM of Al(fcc)100(N x M unit cell)
   4 # =======================================
   5 from Dacapo import Dacapo
   6 from ASE import Atom,ListOfAtoms
   7 from Numeric import sqrt,array
   8 from ASE.Visualization.VMD import VMD
   9 import os
  10 
  11 a0     = 4.05         # cubic fcc lattice constant
  12 N      = 2            # repetition along x
  13 M      = 2            # repetition along y
  14 layers = 2            # slab layers
  15 electronsperatom = 3 
  16 vaclay = 5            # interlayer dist = a0/2
  17 
  18 
  19 # ---------------------------------------------------
  20 # ------- Setup the structure/calculation -----------
  21 # ---------------------------------------------------
  22 
  23 atoms   = ListOfAtoms([])
  24 for n in range(layers):
  25     for i in range(N):
  26         for j in range(M):
  27             scaledpos = [(i+(n%2)/2.)/sqrt(2.),(j+(n%2)/2.)/sqrt(2.),-n/2.]
  28             atoms.append(Atom('Al', a0*array(scaledpos)))
  29                                        
  30                                        
  31 unitcell = [[N/sqrt(2.), 0.0,        0.0],
  32             [0.0,        M/sqrt(2.), 0.0],
  33             [0.0,        0.0,        (vaclay+layers)/2.]]
  34 
  35 atoms.SetUnitCell(a0*array(unitcell),fix=True)
  36 
  37 # VMD(atoms)
  38 
  39 # check if we have already calculated the nc file
  40 if not os.path.isfile('Al100.nc'):
  41 
  42     calc = Dacapo(planewavecutoff = 150, 
  43               nbands = 10 + len(atoms)*electronsperatom/2, 
  44               kpts=(4,4,1),
  45               xc = 'LDA',
  46               usesymm=True,
  47               out = 'Al100.nc',txtout = 'Al100.txt') 
  48 
  49     atoms.SetCalculator(calc) 
  50 
  51 else: 
  52     atoms = Dacapo.ReadAtoms('Al100.nc')
  53 
  54 energy = atoms.GetPotentialEnergy()
  55 
  56 # =======================================
  57 #    STM of Al(fcc)100(N x M unit cell)
  58 # =======================================
  59 # --------------------------------------------------
  60 # ------------ STM image generation part -----------
  61 # --------------------------------------------------
  62 
  63 from Dacapo.ElectronicStates import ElectronicStates
  64 from Dacapo import Dacapo
  65 from ASE.Utilities.STMTool import STMTool
  66 from ASE.Utilities.STMLineScan import STMLineScan
  67 
  68 atoms = Dacapo.ReadAtoms('Al100.nc')
  69 from ASE.Visualization.VMD import VMD
  70 # VMD(atoms)
  71 
  72 loe = ElectronicStates(filename='Al100.nc')
  73 stm = STMTool(loe,contourvalue=1.0e-4, channel="s",
  74              normalaxis=2, smoothfactor=0.1)    # available channels: s,px,py,p
  75 
  76 
  77 
  78 from pylab import *
  79 import matplotlib
  80 
  81 # select a linescan
  82 linescan = STMLineScan(stm,fftnumber=12,axis=1)
  83 
  84 # contour surface plot 
  85 extent = stm.GetExtentOfContourSurface()
  86 im2 = imshow(stm.GetContourSurface().Array,
  87              interpolation='bicubic',
  88              origin='lower',cmap=cm.jet, 
  89              extent = extent) 
  90 colorbar()
  91 
  92 # add position of line scan to contour plot
  93 linex,liney  = linescan.GetLine()
  94 plot(linex,liney,linewidth=2,color='black')
  95 axis(extent)
  96 
  97 savefig('contourplot')
  98 
  99 # line scan plot
 100 figure()
 101 
 102 values = linescan.GetArray()
 103 plot(values)
 104 xlabel('Distance along surface')
 105 ylabel('Height above surface')
 106 title('Linescan at FFT number 12 along the y-axis') 
 107 grid(True)
 108 
 109 savefig('linescan')