Tutorial 3 STM images - Al(100)

The STM is a revolutionary experimental surface probe that has provided direct local insight into the surface electronic structure. Sometimes the interpretation of STM topographs are not straightforward and therefore theoretically modeled STM images may resolve conflicting possibilities and point to an underlying atomistic model. The CAMPOS code includes python modules for generating Tersoff-Hamann STM topographs. Using this functionality requires VTK installed. The STM code is illustrated here for a Al(100) in a 2x2 unit cell.

The setup for the atoms in the Al(100) in a 2x2 unitcell is described below.

After the import statements, the section starting with a0 = ... defines some major variables in the script ; this is not at all required, but increases the readability and reusability of the script. The part starting at atoms=ListOfAtoms(..) illustrates a three-fold python loop that iteratively builds up the slab in a layerwise fashion. After this, atoms.SetUnitCell(...) attaches the periodic boundary conditions to the structur:

>>> from ASE import ListOfAtoms, Atom
>>> from Numeric import sqrt,array
>>> a0     = 4.05         # cubic fcc lattice constant
>>> N      = 2            # repetition along x
>>> M      = 2            # repetition along y
>>> layers = 2            # slab layers
>>> electronsperatom = 3
>>> vaclay = 5            # interlayer dist = a0/2
>>> atoms   = ListOfAtoms([],periodic=True)
>>> for n in range(layers):
>>>     for i in range(N):
>>>         for j in range(M):
>>>             scaledpos = [(i+(n%2)/2.)/sqrt(2.),(j+(n%2)/2.)/sqrt(2.),-n/2.]
>>>             atoms.append(Atom('Al', a0*array(scaledpos)))
>>> unitcell = [[N/sqrt(2.), 0.0,        0.0],
...            [0.0,        M/sqrt(2.), 0.0],
...            [0.0,        0.0,        (vaclay+layers)/2.]]
>>> atoms.SetUnitCell(a0*array(unitcell),fix=True)

Now a calculator must be defined, in this tutorial we will make a STM image from both the GPAW real space grid code and from the dacapo planewave code.

For the GPAW code the calculator for the Al(100) surface can be defined like this:

>>> from gpaw import Calculator
>>> calc = Calculator(gpts=(20,20,48),nbands=28,
...                   kpts=(4,4,1),out='Al100.out'))
>>> atoms.SetCalculator(calc)
>>> energy = atoms.GetPotentialEnergy()
>>> calc.Write('Al100.nc')

The calculator for the dacapo code is similar:

>>> from Dacapo import Dacapo
>>> calc = Dacapo(planewavecutoff = 150,nbands = 28,
...               kpts=(4,4,1),xc = 'LDA',
...               usesymm=True,
...               out = 'Al100.nc',txtout = 'Al100.txt')
>>> atoms.SetCalculator(calc)
>>> energy = atoms.GetPotentialEnergy()

Now a ElectronicStates object must be defined from the output files. This can be done like this for the GPAW code:

>>> from ASE.Utilities.ElectronicStates import ElectronicStates
>>> electronicstates = ElectronicStates(filename = 'Al100.nc')

and for the dacapo code:

>>> from Dacapo.ElectronicStates import ElectronicStates
>>> electronicstates = ElectronicStates(filename='Al100.nc')

dacapo must import its own version of the ElectronicStates class, while the GPAW code can use the generic ASE version.

Now generate the STMTool object from the electronicstates object:

>>> from  ASE.Utilities.STMTool import STMTool
>>> stm = GetSTMTool(electronicstates,contourvalue=1.0e-6, channel="s", normalaxis=2)

The initialization is done with three keyword arguments (keyword=actualvalue) - more options are possible, see section STMTool in the manual.

The most important options are:

The density of states at the Fermi level for which the topograph is generated. This argument is mandatory, as no universal (simple) default is natural. Often used values are [1.0e-6 - 1.0e-3], but is very system dependent and requires attention.
Selects whether the topograph is generated by the value or gradient isosurface of the Fermi level DOS. This adjusts for the orbital character of the DOS of the actual STM tip in the probing energy window, i.e. the STM tip states that couples to the substrate. Available options are "s", "px", "py" and "pz". Default is "s".
The surface normal direction, (0,1,2) for (x,y,z) respectively. normalaxis = 2 is default.
Thermal convolution applied to the slab states. It is necessary to apply this due to finite k-point sampling.

You may additionally initialize the STM tool with voltage for modeling a non-infinitesimal bias voltage on the STM tip, for more details see the section STMTool in the manual.

Finally, the lines:

>>> from ASE.Visualization.VTK import VTKPlotSTM
>>> stmplot = VTKPlotSTM(stm)
>>> stmplot.SetColorRange((3.4,3.8))    # (3.4,3.8) AA distance span color range
>>> stmplot.SaveAsBMP("Al100_stm.bmp")

pops up the STM topograph in a VTK vindow and saves it in bitmap format for you.

>>> stmplot.SetColorRange(((3.4, 3.8))    # (3.1, 3.6) for the dacapo case

will set the color map, so that the electron density at 3.4 Angstroem along the normal axis (starting from the unit cell origin) will map onto black color, whereas the other extreme color, yellow, will be mapped onto the electron density at 3.8 Angstroem along the normal axis. The contour color order is black-red-yellow, so that it resembles the intuitive thermal radiation color scale. Please notice that the STM generation tool is not very automaticed: you need to play around with the STM generation parameters discussed above - STM topograph modeling is an interactive process. Also be aware that Tersoff-Hamann STM pictures for some systems, e.g. gold, are in less good agreement with experimentally obtained STM pictures.

The STM topographs from the GPAW and the dacapo code are shown below:

stm1 stm2

You can combine the STM plot with the atoms, by adding a VTK atoms plot to the stmplot using:

>>> from ASE.Visualization.VTK import VTKPlotAtoms
>>> atoms = Dacapo.ReadAtoms('Al100.nc')
>>> atomplot = VTKPlotAtoms(atoms,parent=stmplot)
>>> stmplot.Render()                        # updates picture

The resulting image is shown below:


Linescans and matplotlib

In this section we will make simulated STM linescans and contour plot using matplotlib.

The example below make a contour plot and a linescan plot. The position of the linescan is plotted on the contour plot. The simulated linescans are made using the STMLineScan class:

import matplotlib
from pylab import *

# select a linescan
linescan = STMLineScan(stm,fftnumber=12,axis=1)

# contour surface plot
extent = stm.GetExtentOfContourSurface()
im2 = imshow(stm.GetContourSurface().Array,
             extent = extent)

# add position of line scan to contour plot
linex,liney  = linescan.GetLine()


# line scan plot

values = linescan.GetArray()
xlabel('Distance along surface')
ylabel('Height above surface')
title('Linescan at FFT number 12 along the y-axis')


The resulting images are shown below:

Contour plot with the position of the linescan added. Simulated STM linescan.