A Guide to simulation of STM images and

spectra from first principles: bSKAN 3.6

W. A. Hofer
Surface Science Research Centre,
The University of Liverpool, Liverpool L69 3BX

August 25, 2005

2
Contents

1 Versions 9
1.1 Additions version 3.6 (2005) . . . . . . . . . . . . . . . . . . . . . 9
1.2 Additions version 3.5 (2004) . . . . . . . . . . . . . . . . . . . . . 9
1.3 Additions version 3.4 (2004) . . . . . . . . . . . . . . . . . . . . . 9
1.4 Additions version 3.3 (2003) . . . . . . . . . . . . . . . . . . . . . 10
1.5 Additions version 3.1 (2003) . . . . . . . . . . . . . . . . . . . . . 10
1.6 Additions version 2.1 (2000) . . . . . . . . . . . . . . . . . . . . . 10
1.7 Original version 1.0 (1999) . . . . . . . . . . . . . . . . . . . . . . 10

2 Introduction 11
2.1 Suitability and systems . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Help utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.1 Keywords and input . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Input errors . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Copyright and license issues . . . . . . . . . . . . . . . . . . . . . 12

3 THEORETICAL BACKGROUND 15
3.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.1 Bardeen approach . . . . . . . . . . . . . . . . . . . . . . 15
3.1.2 Scattering method . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2.1 Topographic mode . . . . . . . . . . . . . . . . . . . . . . 18
3.2.2 Spectroscopies . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Interfaces to DFT programs . . . . . . . . . . . . . . . . . . . . . 21

4 INSTALLATION 23
4.1 Parallel version . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Input features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

5 PROGRAM EXECUTION 25

6 TOPOGRAPHIES 27
6.1 Tersoff-Hamann method . . . . . . . . . . . . . . . . . . . . . . . 27
6.1.1 Possible errors . . . . . . . . . . . . . . . . . . . . . . . . 28

3
4 CONTENTS

6.2 Bardeen method . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6.2.1 Possible errors . . . . . . . . . . . . . . . . . . . . . . . . 29
6.3 Magnetic surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

7 SPECTROSCOPIES 31
7.1 Tersoff-Hamann model . . . . . . . . . . . . . . . . . . . . . . . . 31
7.2 Bardeen method . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.3 Magnetic surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.4 Differential Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . 32

8 EVALUATION 35
8.0.1 Creating current maps and current contours . . . . . . . . 35
8.1 Corrugation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
8.2 Two dimensional maps in parallel planes . . . . . . . . . . . . . . 36
8.3 I/V spectra of a surface . . . . . . . . . . . . . . . . . . . . . . . 36
8.4 Magnetic calculations . . . . . . . . . . . . . . . . . . . . . . . . 36

9 FILES 39
9.1 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
9.2 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

10 KEYWORDS 41
10.1 METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
10.1.1 TERSOFF-HAMANN . . . . . . . . . . . . . . . . . . . . 41
10.1.2 STERSOFF . . . . . . . . . . . . . . . . . . . . . . . . . . 41
10.1.3 NUMERICAL . . . . . . . . . . . . . . . . . . . . . . . . 42
10.1.4 SPECTROSCOPY . . . . . . . . . . . . . . . . . . . . . . 42
10.1.5 FORCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
10.1.6 WAVE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
10.2 SETUP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
10.2.1 ANTIFERROMAGNETIC . . . . . . . . . . . . . . . . . 43
10.2.2 AREA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
10.2.3 BIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
10.2.4 CELL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
10.2.5 DELTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
10.2.6 FERROMAGNETIC . . . . . . . . . . . . . . . . . . . . . 45
10.2.7 GRIDPOINTS . . . . . . . . . . . . . . . . . . . . . . . . 45
10.2.8 HOLLOW . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
10.2.9 LIMITS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
10.2.10 NKELDYSH . . . . . . . . . . . . . . . . . . . . . . . . . 46
10.2.11 NSPECTRUM . . . . . . . . . . . . . . . . . . . . . . . . 46
10.2.12 PIVOT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.2.13 TOP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.2.14 ZVACUUM . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.3 EVALUATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.3.1 CURRENT . . . . . . . . . . . . . . . . . . . . . . . . . . 48
CONTENTS 5

10.3.2 CORRUGATION . . . . . . . . . . . . . . . . . . . . . . . 49
10.3.3 MERGE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
10.3.4 PLOTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

11 Spectroscopy evaluations 51

12 WAVEFUNCTIONS 53

13 Geometry files 55
6 CONTENTS
Preface

This guide is intended as a hands on manual for the execution of bSKAN 3.6, an
optimized and parallel code to simulate STM topographies and spectroscopies
from first principles. The program is an open source package, it can be used
free of charge. However, use of the program is limited to users complying with
two conditions: (i) Acknowledgement of the source; and (ii) feeding back all
improvements made to the code to the original author. These conditions are
mandatory, and users who are found not to comply with these rules will be
excluded from future releases. The program requires a minimum of 2GHz pro-
cessors, with a memory of no less than 1GB. In parallel mode it has been tested
for up to 64 processors. The memory requirement for high level computations
of systems of medium size is about 200-500 MB.

SKAN b 7
8 CONTENTS
Chapter 1

Versions

1.1 Additions version 3.6 (2005)

The main changes were in (i) the rewrite of the program to account for the newly
developed theoretical method based on the Keldysh formalism, (ii) changes in
the evaluation routines, and (iii) rewriting the symmetry analysis and the gener-
ation of symmetry, working now automatically. The default calculation is now
with the standard Bardeen method (NKELDYSH = -1), the bias dependent
corrections are computed with NKELDYSH = 1. The program now integrates
the differential maps at the end of the calculation, the evaluation then can be
performed with a setpoint taken from experiments (I,V values).

1.2 Additions version 3.5 (2004)

The main changes are the spectroscopy modules of bSKAN 3.5. Model cal-
culations showed that (i) the reduced number of layers in the tip description,
and (ii) the numerical stability on parameters for spectroscopy calculations is
not satisfactory, if spectra are obtained by a numerical differentiation of I(V)
maps. Version 3.5 therefore performs spectroscopy calculations differentially,
the new routines allow to identify unambiguously the effect of the tip electronic
structure.

1.3 Additions version 3.4 (2004)

Additional module for differential spectroscopy, in version 3.4 only as addition to
the numerical differentiation of I(V) maps. Spectroscopy module optimized for
calculations with more than thousand k-points per system. Included a method
to calculate chemical interactions and their effect on tunneling currents, as de-
scribed in PRL 92, 266101 (2004).

9
10 CHAPTER 1. VERSIONS

1.4 Additions version 3.3 (2003)

Much improved version of data representation by creating an interface which
can be visualized with OpenDx throughout all calculation methods.

1.5 Additions version 3.1 (2003)

Tunneling topography and spectroscopy now equally implemented. Massive
code optimization increased the speed of calculations by two orders of magni-
tude, which is particularly important for complicated metal systems. Parallel
code optimized.

1.6 Additions version 2.1 (2000)

Optimization of the code and first high resolutions simulations. Parallel code
developed.

1.7 Original version 1.0 (1999)

Based on FLAPW wavefunctions and only serial, this very first version could
only calculate a few points on metal surfaces, with a rather limited resolution
of k-space.
Chapter 2

Introduction

The program bSKAN is written in modular form and in the current imple-
mentation in Fortran 90. In contrast to Fortran 77 this allows to use derived
memory structures like types, which in turn make the allocation of memory, the
transfer of data, and the handling of large and complex structures much easier.
For example, the whole package is programmed without a single common block,
and all sizeable memory is allocated during runtime.
The main programming challenge was the reduction of operations. As the
wavefunctions are given in a two dimensional Fourier grid of typically more than
one hundred components, the integration involves handling of matrices of ten
thousand components. This can only be accomplished in a reasonable timescale
if all steps are highly optimized.
At present, and in a serial implementation, the program is a able to calculate
a single gridpoint of the STM tip position in timescales of typically less than
one minute, which makes the calculation of detailed images in high resolution
possible within a few hours. In parallel execution e.g. on a SGI R10000 cluster,
we have calculated the spectrum of one point on a magnetic surface, mapped
with 3000 k-points in the IBZ, using a model tip with 400 k-points around the
centre of the Brillouin zone, and an energy grid from -1V to +1V of 101 points
in less than four hours. The program in this case uses about 2GB of memory.

2.1 Suitability and systems

So far simulations of STM and STS have been performed on a wide variety of
systems: magnetic and non-magnetic metals, semiconductors, semiconductors
with magnetic properties, molecules on metal and semiconductor surfaces, oxy-
gen covered metals etc. In all cases the simulations agree reasonably well (same
order of magnitude for the current in simulations and experiments for a given
result) to spectacularly well (same current values). The qualifying facts for a
given calculation seem to be: (i) Whether all effects are included in groundstate
DFT calculations (here one can be sceptical, in particular if highly correlating

11
12 CHAPTER 2. INTRODUCTION

systems are analyzed), and (ii) whether the experimental range is reasonable for
perturbation methods (here, as a rule of thumb, we are limited to a maximum
of about 5-10nA on metals for low voltages). Within this range the calculations
should be generally safe and easy to perform.

2.2 Help utilities

2.2.1 Keywords and input
From the viewpoint of users it seemed important to structure the input in an
easy manner. The main input is therefore reduced to a limited number of
keywords (see Appendix), and the input routine provides help functionalities
for input errors. The easiest way to get started is to provide an input file with a
single line HELP, which invokes a routine writing a file README detailing
all the options.

2.2.2 Input errors

The program contains a rudimentary - far from complete - check of input data
for plausibility. In every run, where a problem is detected, a README file is
created which specifies the problem. Mainly these are:

• The energy range of eigenvalues is smaller than the energy range of the
calculation. Remedy: go back to the DFT calculation where the STM
wavefunctions were created and increase the energy interval.

• The energy resolution of a spectroscopy calculation is too high for the

input wavefunctions. This is usually correct for metals, where the eigen-
values are densely spaced, but not necessarily correct for semiconductors,
where you have a bandgap. The routine checks whether every interval con-
tains at least one eigenvalue. Remedy: increase the number of k-points.

• The k-points of either surface or tip do not cover the IBZ of the first
Brillouin zone. This can be intended, if for example only a limited region
of the IBZ is considered, or it can be an error, if the k-point sampling is
incomplete. Remedy: check and if necessary change the k-point sampling
in the DFT calculation.

2.3 Copyright and license issues

The copyright of the program rests with the authors. However, the program
is distributed as open source program. This means that no licence fees apply,
but also, that extensions and improvements of the programs should be made
available to other users. Please include a reference to [1, 2] in every work which
uses the Bardeen integration, a reference to [3], if you calculate spin-resolved
currents, and a reference to [4], if you perform spectroscopy simulations. The
2.3. COPYRIGHT AND LICENSE ISSUES 13

extension to multiple scattering in the vacuum barrier is based on a publication

in Journal of Physics [5].
14 CHAPTER 2. INTRODUCTION
Chapter 3

THEORETICAL
BACKGROUND

STM theory, like STM experiment, has a history of at least twenty years, from
the earliest papers of Binnig and Rohrer [6, 7, 8] on the origin of the instru-
ments precision to the theoretical models of Büttiker and Landauer [9], Tersoff
and Hamann [10, 11], Chen [12, 13], Sautet and Joachim [14], or Flores [15, 16].
Basically, there are two different philosophies concerning the importance of dif-
ferent effects on STM images and spectra:
On the one hand, it was thought that the scattering process itself contains
the main physical parameters determining the images. This is reflected in all
scattering approaches, where the exact electronic structure of the two surfaces
is commonly treated in a very rudimentary fashion.
On the other hand, it is thought that the scattering process, due to the
large timescales involved (the interval from one tunneling electron to the next is
typically larger than picoseconds), makes a scattering approach redundant, as
long as the electronic structures of the two surfaces are well described. bSKAN
follows the second line of argument. The wavefunctions of both surfaces are
determined by highly precise density functional calculations, while the transi-
tion process is described by perturbation theory. The theoretical model goes
back to Bardeen’s treatment of a metal-insulator-metal junction [17], and it
has been used in the last years for a wide range of materials from metals and
magnetic overlayers to semiconductors and molecules adsorbed on metals and
semiconductors. For a review see [2].

3.1 Method
3.1.1 Bardeen approach
Within the Bardeen approach to tunneling the current between a surface and a
tip is described by the sum over surface and tip states as follows:

15
16 CHAPTER 3. THEORETICAL BACKGROUND

¯ Z ¯2
4πe X ¯ h̄2
¯ ∗ ∗ ¯
¯
I= [f (Eµ ) − f (Eν + eV )] ¯− dS (χν ∇ψµ − ψµ ∇χν )¯
h̄ µ,ν 2m S
δ (Eν − Eµ + eV ) (3.1)

Here, ψ is the wavefunction of the single electron state (in DFT Kohn-
Sham state) of the surface, χ a single electron state of the tip, f is the Fermi
distribution function, the bias voltage between surface and tip equals V , and
the integration surface is assumed to be in the vacuum region. The key variable
in this relation is the integral over the separation surface, which is called the
tunneling matrix element Mµν . It is defined by:
Z
Mµν = dS (χ∗ν ∇ψµ − ψµ ∇χ∗ν ) (3.2)
S
The matrix element is a scalar quantity, which is equivalent to the overlap
of the vacuum wavefunctions of surface and tip. To implement this approach
within the periodic systems typical for groundstate DFT calculations, the fol-
lowing points have to be considered:

• The wavefunctions in DFT are given as Kohn-Sham states of specific

points of the two dimensional Brillouin zone of the surface. Each k-point
has its own range of eigenvalues and states.
• Since the lattice geometry of surface and tip are in general incommen-
surate, each k-point of the surface is inequivalent to each k-point of the
tip.
• Most DFT codes reduce the number of operations to achieve convergence
by utilizing symmetry properties of their systems. The k-points of a given
mesh reflect these properties.

The method used in bSKAN accounts for these properties of groundstate

DFT calculations in the following way: (i) The lateral k-value of a given state is
not conserved in the transition. This means that all transitions are admissible
as long as the electron energy is conserved. (ii) The wavefunctions of the DFT
input are expanded over the whole Brillouin zone using the symmetry operations
of the underlying lattice. For a lattice of hexagonal symmetry this means, for
example, that every wavefunction read in is equivalent to six wavefunctions
determined by the rotation of reciprocal lattice vectors.
The wavefunctions required by bSKAN have the following form:
X
ψµ (kk , rk , z) = Cµ (Gk , z) exp i(kk + Gk )rk (3.3)
Gk

At present the grid in z-direction is hardwired in the program at 0.1 a.u.

(0.05218 Å). It was found that this resolution is sufficient to reproduce corruga-
tion values of metal surfaces precisely down to a corrugation amplitude of less
3.1. METHOD 17

than 1 pm, which is about the resolution of todays best instruments. Given
the usual size of systems in DFT (topographies on metals: 10-40 k-points in
the IBZ, six to eight symmetry operations, expansion up to 200 G-vectors), a
point by point integration over the separation surface is ruled out for practical
reasons. Therefore an additional assumption is made in bSKAN: The separa-
tion surface is a plane located in the middle between the tip and the surface.
In this case the integration for a single Fourier component can be performed
analytically, provided the region of the surface is limited. This is generally the
case, if the tip consists not of a plane surface, but a surface with an attached
microtip of one to a few atoms and one or more layers. It is established opinion
today that such a tip is used in all high resolution scans. For calculated model
tips under these assumptions see [2].

3.1.2 Scattering method

From a theoretical point of view a tunneling electron, e.g. in a scanning tun-
neling microscopy measurement, is part of a system comprising two infinite
metal leads and an interface, consisting of a vacuum barrier and, optionally, a
molecule or a cluster of atoms with different properties than the infinite leads.
The system can be said to be open - the number of charge carriers is not con-
stant - and out of equilibrium - the applied potential and charge transport itself
introduce polarizations and excitations within the system. The theoretical de-
scription of such a system has advanced significantly over the last years, to date
the most comprehensive description is based either on a self-consistent solution
of the Lippman-Schwinger equation or on the non-equilibrium Green’s function
approach. Within the vacuum barrier itself, inelastic effects play an insignifi-
cant role. Here, as in most experiments in scanning tunneling microscopy, the
problem can be reduced to the description of the tunneling current between two
leads - the surface S and the tip T - thought to be in thermal equilibrium. The
bias potential of the circuit is in this case described by a modification of the
chemical potentials of surface and tip system, symbolized by µS and µT . This
reduces the tunneling problem to the Landauer-Büttiker formulation:
Z +∞
2e £ ¤
I= dE [f (µS , E) − f (µT , E)] × T r ΓT (E)GR (E)ΓS (E)GA (E)
h −∞

Here, f denotes the Fermi distribution function, GR(A) (E) is the retarded
(advanced) Green’s function of the barrier, and ΓS , ΓT are the surface and
tip contacts, respectively. They correspond to the difference of retarded and
advanced self energy terms of surface and tip; we define them by their relation
to the spectral function AS(T ) of the surface (tip):
h i
AS(T ) (E) = i GR A R A
S(T ) (E) − GS(T ) (E) = GS(T ) (E)ΓS(T ) (E)GS(T ) (E) (3.4)

The multiple scattering formalism can be evaluated in real space, with the
18 CHAPTER 3. THEORETICAL BACKGROUND

help of an eigenvector expansion of the surface and tip Green’s functions:

R(A)
X ψk (r1 )ψ ∗ (r2 )
k
GS (r1 , r2 , E) = (3.5)
E − Ek + (−)iη
k

R(A)
X χi (r1 )χ∗ (r2 )
i
GT (r1 , r2 , E) = (3.6)
i
E − E i + (−)i²

The zero order current results

· µ ¶ µ ¶¸ ¯µ ¶ ¯2
4πe X eV eV ¯
¯ h̄2 eV ¯
¯
I(0) = f µS , Ek − − f µT , Ei + ¯ − 2m − κ2 − κ2 Mik ¯ δ(Ei −Ek +eV ).
h̄ 2 2 i k
ik
(3.7)
The terms κ denote the vacuum decay of the surface (k) and tip (i) wavefunc-
tions.
The result for the first order current, including only the terms for single electron
paths (essentially the square of the matrix Mik , while multiple scattering path-
ways will be described by four and six matrix multiplications), then involves
also a term which depends on the bias voltage:
· µ ¶ µ ¶¸ ¯µ ¶ ¯2
4πe X eV eV ¯
¯ h̄2 eV ¯
I(1) = f µS , Ek − − f µT , Ei + ¯ − + 2 2 Mik ¯¯ δ(Ei −Ek +eV ).
h̄ 2 2 2m κi − κk
ik
(3.8)
It can be seen that the zero and first order currents differ only in the sign of
the explicit bias dependent part. Moreover, the obtained tunneling currents for
higher voltages will increase more than linearly with the applied bias voltage
and for both the special case of zero bias results us exactly the Bardeen current.
The method is presented in [5].

3.2 Implementation
3.2.1 Topographic mode
The method is implemented in the program in the following way: first the lattice
parameters of surface and tip are read in and the two lattices are expanded over
the full Brillouin zone. Then a rectangular grid covering the surface unit cell is
set up. The matrix elements MG,G0 are calculated by analytically integrating
the plane wave components over the surface of the tip unit cell:
Z
MGS GT = dS exp i(GS − GT )r GS/T = kk + Gk (3.9)
S

For a given position R of the STM tip the matrix elements are multiplied
by the phase of the surface wavefunctions:

NGS = exp iGS R (3.10)

3.2. IMPLEMENTATION 19

The current of a given transition µ → ν and at a certain distance d is

therefore the sum over all Fourier components of surface and tip wavefunctions.
It contains three distinct components: (i) The z-dependent amplitudes; (ii) the
integrals and phases depending on the lateral position of the STM tip; and (iii)
the occupation numbers of electrons and a Gaussian, which mimicks the delta
functional of elastic transport, and which depends on the tunneling conditions
and the energy eigenvalues.

4πe X · dCµ (GS , z) dC ∗ (GT , d − z)

¸
Iµν (d) = wµν C ∗ (GT , d − z) − Cµ (GS , z) ν ×
h̄ dz dz z=d/2
GS ,GT
µ ¶
2 (Eµ − Eν + eV )2
× |MGS ,GT NGS | × [f (Eµ − f (Eν + eV )] exp − (3.11)
2σ 2

To speed up the program these three components are calculated separately in

the simulation of topographies. The integral over the tip unit cell of the Fourier
components is calculated initially and stored outside the loop changing the tip
position. The calculation of the energy dependent components is also outside
the loop over the tip positions in topographies. The phases are calculated after
every shift of position of the STM tip. The sum over the z-dependent amplitudes
and derivatives of the wavefunctions is in every case the innermost loop of the
calculation. Simulations are routinely done over the whole range of z-values.
The weight of a given transition wµν depends on (i) the weight of surface and
tip states, and (ii) the decay constants of the surface and tip states. These decay
constants are calculated after the wavefunctions are read and stored in separate
tables, which are used after the integration of Fourier coefficients to determine
the weight of an individual transition. At present the program allows to calculate
the current (topographies) and the differential spectrum (spectroscopies) either
with the standard Bardeen method (no bias dependency), or with zero order or
first order scattering methods, see section Method in this chapter.

3.2.2 Spectroscopies
Initially, spectroscopy functionality was built into the code by an additional
loop over bias voltage. In this case the obtained results were topographies for
every given bias voltage in a bias interval, e.g. from -1V to +1V. A comparison
with experimental spectra (dI/DV spectra) was then simulated by a numerical
differentiation of I(V) for every single point of the surface. The method proved
to have several methodical problems:

1. The change of occupation numbers near the Fermi level due to the Fermi
distribution function always shows up as a distinct spike in the spectrum.

2. The number of layers of the surface and tip electronic structure deter-
mines the spacing of eigenvalues and thus the ensuing spectrum (minimum
number for noble metals about 23 layers). This is close to impossible to
20 CHAPTER 3. THEORETICAL BACKGROUND

calculate for the tip electronic structure, because the tip requires a very
large unit cell of at least eight atoms per plane.
3. The ensuing I(V) curve and its numerical derivative do not allow a clear
identification of surface and tip contributions and make it very difficult,
if differences to experimental spectra are observed, to improve the repre-
sentation.

For these methodical reasons the frontal attack of the problem was finally
given up, after a large number of trial calculations on Fe, Cr, Mn/Fe systems
and Cu, Ag, and Au surface states. Instead, the program now contains a com-
prehensively differential approach to the problem. The details of the theoretical
analysis and the new approach are published in [4]. For the present purpose
the relevant result is that the differential spectrum dI/dV is directly calculated
and written to a file, which contains the dI/dV map for a defined surface grid.
From this spectrum the I(V) map is obtained by integration. The incremental
change in the current due to a change of bias from V to V + dV is:
X X
dI = |M (ψi1 , χk1 )| + |M (ψi2 , χk2 )| (3.12)
i1 k1 i2 k 2

where the eigenvalues of surface Ei1(2) and tip Ek1(2) states are within the
intervals:

Ei1 ∈ [EF + eV − edV /2, EF + eV + edV /2]

Ek1 ∈ [EF − edV /2, EF + edV /2]
Ei2 ∈ [EF − edV /2, EF + edV /2]
Ek2 ∈ [EF − eV − edV /2, EF − eV + edV /2] (3.13)

Here, EF denotes the Fermi level of surface and tip system, respectively.
Then the total spectrum contains equally two distinct contributions due to the
bandstructure of the surface and the tip system:

dI(V ) X |M (ψi1 , χk1 )| X |M (ψi2 , χk2 )|

= + (3.14)
dV dV dV
i1 k1 i2 k 2

The files containing the two separate contributions to the spectrum and the
sum of these contributions. It is therefore possible to identify the origin of a
feature in the spectrum and determine, whether it is due to the surface or the tip
electronic structure. In general we find that surfaces with a very low density of
states at the Fermi level lead to spectra with features unique of the surface. This
is for example true for semiconductor surfaces. For metal surfaces we find that
this is not the case and the ensuing spectra will show up both bandstructures. In
this case there are two options to obtain a reliable result: (i) Limit the analysis
to the surface contribution alone, which maps to the states at the Fermi level
of the tip; or (ii) or increase the number of layers and the precision of the tip
3.3. INTERFACES TO DFT PROGRAMS 21

model so that its bandstructure is correctly represented in the calculation. In

principle, this is feasible and will become routine once computing power has
increased to this level. For the time being, we suggest to use the first approach.

3.3 Interfaces to DFT programs

The Kohn-Sham states can in principle be obtained from any DFT method,
which describes the electron states in the vacuum range by a two dimensional
Fourier expansion (x,y) and a real space grid (z). At present interfaces ex-
ist for VASP (see [Link] and film FLAPW. Since
FLAPW has been substantially altered during the last years by Michael Weinert
and Raimund Podloucky, its interface is no longer up to date and will have to be
rewritten in the future. It is also possible to use other methods like SIESTA or
Wien2k for this purpose, specifications for the interfaces can be obtained from
WH.
22 CHAPTER 3. THEORETICAL BACKGROUND
Chapter 4

INSTALLATION

The program is delivered in a compressed [Link] file, which needs to be uncom-

pressed either with gunzip (the Unix utility), or any other of available utilities
like WinZip etc. Once the Fortran 90 modules (named *.F) and the make-
files are stored in a directory, the executable can be compiled with any suitable
compiler. The current implementation supports Intel Fortran Compilers (Make-
[Link]), Sun clusters ([Link]) and Silicon Graphics clusters ([Link]).
Portation to other systems should be unproblematic, because the program does
not depend on external libraries. All routines use standard Fortran 90. The
commands to unzip and to extract the program files are:
gunzip [Link] tar -vxf [Link]
After this the [Link] needs to be copied onto makefile, e.g. for a Silicon
Graphics environment with:
cp [Link] parallel makefile
Then the (parallel) executable can be compiled with:
make bskan
For serial executables other makefiles with the tag serial have to be used.
Please note that the location of the libraries as well as the switches generally
depend on the setup of your cluster. The easiest way to find out about your
environment is to ask the system administrator. There should be no error
messages during compilation. If, however, the compilation ends with an error,
please check first that the correct makefile was used. If this is the case and the
errors do not disappear, please contact your system administrator.

4.1 Parallel version

The problems with parallel coding are well known: there exists no standard
implementation of the MPI interface e.g. for Linux clusters. The necessary
libraries (LAPACK, BLAS, SCALAPACK) need to be compiled for the com-
puter environment. However, this is not usually trivial and will best be done by
the administrator. I therefore follow the usual practice not to provide explicit

23
24 CHAPTER 4. INSTALLATION

advice for MPI implementation. The parallel compilation, as provided for the
computer environment the code has been running so far, can be seen in the
makefiles. The location of the libraries will invariably vary, as will the compiler
used (for example, I used a PGF90 compiler on Linux 2.7.2 clusters, built from
AMD Athlon processors). The makefile therefore has to be modified. Again,
this is best done in cooperation with the computer administrators.

4.2 Input features

Generally, it was sought to minimize the input to the bare essentials for a run.
The program therefore provides a number of default settings, which are written
at the beginning of the output file. It was also sought to make the input format
as free as possible. However, there is a limit, where coding becomes rather
demanding, without a substantial gain in efficiency or flexibility. The input
routines of bSKAN are aimed at a compromise: a command usually consists of
one word plus zero or more values. The first command in the input file should
always be the method command. Apart from that the order of commands is
optional.
Chapter 5

PROGRAM EXECUTION

The name of the executable is bskan36, it can be either executed in interactive

mode by bskan36 or bskan36 &. Initially, the program searches in the same di-
rectory for five files: INSCAN, WAVSAMPLE, WAVTIP, ASAMPLE and ATIP.
The first contains all the input parameters, the second two the wavefunctions
of the sample surface and tip, respectively, and the last two contain the atomic
positions in direct coordinates for the surface and tip, respectively (see section
FILES). Please note that WAVTIP and ATIP files are not needed in case of a
Tersoff-Hamann calculation.
For most applications the program will be executed in a queueing system.
Please remind that it might take some time (up to a few minutes, depending on
the system), for the program to produce any output. In case the program stops
by writing a README file, it detected some input error in the file INSCAN.
In case it stops without such a message, the file OUTSCAN should contain a
message detailing an error in a read operation on the wavefunction files. In this
case the files are probably corrupt and have to be generated again. Depending
on the tasks defined in INSCAN, the program generates various output files.
These are listed in the section FILES.

25
26 CHAPTER 5. PROGRAM EXECUTION
Chapter 6

TOPOGRAPHIES

In topographic mode an STM scans across the surface while the tunneling cur-
rent is kept constant. This could be mimicked by a suitable feedback within the
program. However, the feedback algorithm is tricky to program and convergence
then becomes a major issue. Therefore it was decided to compute a complete 3D
matrix of tunneling currents on the surface. This is still manageable, computa-
tionally, and it makes the program substantially simpler. The vertical extension
depends on the input files, generally all z values of the smaller file are included.
If, for example, the WAVSAMPLE file includes 50 z-values, and the WAVTIP
file 100, then only the first 50 gridpoints of the tip electronic structure are
included in the calculation. bSKAN provides two different routines:

6.1 Tersoff-Hamann method

Since this is included in practically all DFT codes, it also provides a handy check
of the surface electronic structure and subsequent calculations with a model tip.
The input in the file INSCAN is the following:
TERSOFF HAMANN MODEL
BIAS VOLTAGE = -0.01
LIMITS = -0.05 0.05
GRIDPOINTS = 61
CELL = 1.0 1.0
PIVOT POINT = 0.0 0.0
NKELDYSH = 1
ZVACUUM = 11.2
BIAS VOLTAGE and LIMITS values are in eV. The limit describes an ambi-
ent environment, it determines the states included in the summation outside the
energy window defined by the bias voltage. The number of gridpoints applies to
the longest axis. For a square lattice, this means a quadratic grid of the surface
mesh. For a rectangular lattice the shorter direction is covered by proportion-
ally less gridpoints, so that the mesh is equally spaced in both directions. For

27
28 CHAPTER 6. TOPOGRAPHIES

a hexagonal lattice it creates a rectangular mesh, where the rectangular lattice

vectors (ARs) are
AR1 = A1 + A2
AR2 = A1 - A2,
with A1 and A2 being the lattice vectors of the hexagonal lattice. It is
easily seen that the rectangular cell has an area of two hexagonal unit cells.
The default without an input of GRIDPOINTS is 31. The image size can be
varied with the keyword CELL. The bSKAN default here is one rectangular
unit cell. The keyword PIVOT POINT determines the lower left point of the
created surface image. The default, if no value is given, is the point (-0.5,-0.5).
NKELDYSH determines whether the scattering approach is used. Finally, the
ZVACUUM parameter describes the vacuum boundary of the sample surface (in
Å), in the above example it is 11.2 Å. For more details on keywords, see chapter
KEYWORDS. The program creates two output files: OUTSCAN provides the
information about the system and the run, the file CURMAT contains the binary
3D matrix of local density of states. The file CURMAT can be used to evaluate
the relevant properties like surface CORRUGATION, the apparent height of
atoms on this surface, it can also be used to create constant density or constant
height contours which can be compared to the experiments. This is described
in the chapter EVALUATION.

6.1.1 Possible errors

There are essentially two classes of errors: either the program stops, without
creating the matrix, or it creates the matrix but gives unexpected results. In
the first case it will either create a README file, then the parameter input
contained an error and the file should contain information to correct the error.
Or if it does not, then the wavefunction file is corrupted and should be generated
again. In the second case the possible sources of error are the input range of
eigenvalues. If the limits and the bias voltage lead to an energy range which
is beyond the limits defined in the wavefunction file, the result will be an error
message but no stop of the program.

6.2 Bardeen method

In this method the tip is included in the calculation. Two files, WAVTIP and
ATIP are needed in the working directory containing the electronic structure of
the STM tip model, and the atomic positions of the tip in direct coordinates,
respectively. The input parameters are the following:
NUMERICAL EVALUATION
BIAS VOLTAGE = -0.01
LIMITS = -0.05 0.05
GRIDPOINTS = 61
CELL = 1.0 1.0
PIVOT POINT = 0.0 0.0
6.3. MAGNETIC SURFACES 29

NKELDYSH = 1
ZVACUUM = 11.2
The program in this case sets up a surface grid, computes the matrix of
integrated Fourier components, and determines the eigenvalues to be included
before looping over the surface grid. The time needed for a gridpoint scales
nearly linearly with the bias voltage, since this determines the number of in-
cluded states. For low bias scans on metal surfaces it is commonly in the range
of seconds, for semiconductors the duration is considerably higher and can be as
long as a few minutes. The program produces the usual output files OUTSCAN
and CURMAT, and in addition a formatted file CURSAVE, which should be
saved, since binary files like CURMAT do not always port easily from one sys-
tem to the other. By playing around with different model tips it can be seen
that the Bardeen integration makes tunneling topographies tip dependent to
quite a high degree. This means that the inclusion of the tip adds an additional
dimension in the comparison between experiments and simulations. A good
agreement between them requires that the experimental input (current, bias
voltage), and output (corrugation, shape of a structure) agrees with the input
and output in the simulation.

6.2.1 Possible errors

The energy range of tip and surface electronic structures determines the range
of possible bias potentials. A bias potential outside the range of eigenvalues of
either surface or tip will lead to an error message in the output.

6.3 Magnetic surfaces

While tunneling currents into paramagnetic tips made of tungsten are not spin-
selective - both electron states of a magnetic surface tunnel into the same states
of the tip -, the situation changes for magnetic tips. Here, the spin-up and
spin-down states find a different electronic structure, with commonly a higher
density of spin-down states at the Fermi level. This favors transitions of spin-
down electrons, which leads to a magnetic image of the surface, or an image,
predominantly, of the electronic structure of the minority band. bSKAN in-
cludes functionalities both, for the calculation of spin-polarized currents, and
for the evaluation of contours if the magnetization direction of surface and tip
are not collinear. To make a non-collinear calculation it is first necessary to
determine the currents for both ferromagnetic (up states into up states) and
antiferromagnetic transitions (up states into down states). This is done by
adding the keyword:
FERROMAGNETIC (ANTIFERROMAGNETIC) ORDERING
The two runs will yield two different current matrices, which are merged, in a
second step, under the assumption of an angle between the magnetization vector
of surface and tip (see [2]). This step is described in the section EVALUATION.
30 CHAPTER 6. TOPOGRAPHIES
Chapter 7

SPECTROSCOPIES

In the spectroscopic mode the STM tip is stabilized at a point above the surface,
this point is usually described by a bias voltage/current combination. After sta-
bilization the feedback loop is disengaged, and the bias voltage ramped from a
lower limit to an upper limit. The current/voltage curves in this case look rather
bland, but their first and second derivatives contain information about the sur-
face electronic structure (e.g. surface states on (111) noble metals or (100)
iron) and the dynamic changes due to electron-electron and electron-phonon in-
teractions. Within bSKAN a spectroscopy is a topography over a range of bias
voltages. This makes it possible to study not only the spectroscopies at fixed
points of the surface, but also to study their change with the STM tip position,
which in turn can yield valuable information about local electronic properties.

7.1 Tersoff-Hamann model

The input is similar to the input used for topographies and with the TH-method.
The minimum input for a spectrum is the following:
STERSOFF = -1.0 1.0
LIMITS = -0.05 0.05
NSPECTRUM = 101
GRID = 1
NKELDYSH = 1
ZVACUUM = 11.2
Here, the spectrum covers the interval from -1.0 to 1.0 Volt, the surface
is probed at only one gridpoint (the TOP point, which is (0.0,0.0) in default,
and the energy interval from -1.0 to +1.0 eV is computed with 101 values.
The variables NSPECTRUM and GRID could in principle also be omitted, the
defaults within bSKAN are 11 energy gridpoints (NSPECTRUM) and 31 surface
gridpoints (GRID) along the major axis. NKELDYSH determines whether the
scattering approach is used, since ZVACUUM sets the vacuum boundary. For
more details on keywords, see chapter KEYWORDS. The output of such a run

31
32 CHAPTER 7. SPECTROSCOPIES

consists of three files. OUTSCAN gives, as usual, the information about the
system and the tunneling parameters as well as the included states. The files
CURSPEC and CURSAVE contain the current matrix for all local and energy
gridpoints, the faster loop in this case runs over the energies.

7.2 Bardeen method

The only difference is the method keyword, which has to be changed. The input
for a spectroscopy calculation with the Bardeen method is the following:
SPECTROSCOPY = -1.0 1.0
LIMITS = -0.05 0.05
NSPECTRUM = 101
GRID = 1
NKELDYSH = 1
ZVACUUM = 11.2
The output is the same as above. It is recommended to save the file CUR-
SAVE since the binary CURMAT file is not generally transferable to other
platforms.
In case of spectroscopies the representation of the bandstructure becomes the
most important parameter for the quality of the spectrum. In general, a too low
number of k-points leads to a loss of resolution and even to a complete distortion
of the spectrum. The necessary number of k-points depends to some extent on
the desired resolution, i.e. the energy grid in the calculation. To analyze the
grid the information about the number of states in every interval are printed
out in the file TRANSLOG. In case the grid is too small, the number of states in
an interval approaches one or even reaches zero. In this case a warning message
is printed in the OUTSCAN file. It is recommended to increase the k-point
sampling until this warning disappears.

7.3 Magnetic surfaces

The only additional information needed is the magnetic ordering. The bSKAN
default is ferromagnetic, the explicit keyword for ferromagnetic and antiferro-
magnetic ordering are the following:
FERROMAGNETIC (ANTIFERROMAGNETIC) ORDERING

7.4 Differential Spectroscopy

In general it is desirable to have a clear representation of the STM tip, which
can be inferred from experimental spectra and gives the correct contributions
to a spectrum over a limited voltage around the Fermi level. However, such
a tip cannot be calculated today even with high performance computers. The
reason is twofold: (i) The electronic bandstructure of a metal film is discrete
due to the limited number of layers in the film. The spacing of the eigenvalues
7.4. DIFFERENTIAL SPECTROSCOPY 33

at a given point of the Brillouin zone reflects this limitation. (ii) The number
of two dimensional k-points is also limited, which reduces the precision of the
bandstructure map also in two dimensions.
34 CHAPTER 7. SPECTROSCOPIES
Chapter 8

EVALUATION

8.0.1 Creating current maps and current contours

The keyword for evaluating the current maps is CURRENT. In connection with
a real number it has two different meanings:
CURRENT = 0.0
will create a file CURRENT, which contains the current map in a format com-
patible with OpenDx (this means that the current is given in nested loop of
three indices).
If the keyword is used with a positive value, e.g.
CURRENT = 1.5
then the program constructs a current contour of the surface with this input
value (generally in nA, for TH topographies in units of the LDOS). This contour
is written to the file PLOTCON, which contains in the first line the information
about the contour maximum and minimum as well as the current value.

8.1 Corrugation
The main information contained in an STM image is the corrugation height, or
the difference of the vertical position of the STM tip between a hollow site and
an on top site. This information can be extracted from the current matrix with
the commands:
CORRUGATION
TOP = 0.0 0.0
HOLLOW = 0.5 0.5
The first keyword defines the task, the other keywords define the position of
the ion and the hollow site. The file written is called PLOTCOR, it contains
the z-dependent current values at both sites, the apparent barrier height due to
the current decay, and the corrugation value in Å.

35
36 CHAPTER 8. EVALUATION

8.2 Two dimensional maps in parallel planes

The current can also be plotted in 2dim maps at preset z values from the surface.
This is done with the following commands:
PLOTS = 5
FPLOT = PLT
ZPLOT = 10 20 30 40 50
The program creates five output files, called PLT.001 to PLT.005, containing
the currents in the parallel planes for z = 10 to z = 50. The first line of each
file gives the distance and the maximum current value in this plane.

8.3 I/V spectra of a surface

The spectra of the surface are contained in the file CURSPEC, they cover all
current values over the surface grid for bias voltages within the predefined range.
To extract the currents for a given position of the tip, the bias voltage and the
current value, at which the tip was stabilized, have to be defined. The command
lines to this end are the following:
BIAS VOLTAGE = - 0.3
CURRENT = 1.5
TOP = 0.0 0.0
This creates four output files. The file PLOTSPC contains the current, the
normalized derivative and the second normalized derivative of the current at
the point TOP. The vertical position of the tip in this case is preset to the value
defined by the bias/current values. In addition, a two dimensional map of cur-
rents, first and second derivatives is written to the files PLOT.01 to PLOT.03,
where the map is determined by the position TOP and the bias/current values
from the input.

8.4 Magnetic calculations

Here, two separate outputs can be created with the keywords FERROMAG-
NETIC or ANTIFERROMAGNETIC. These result of the calculations have to
be merged under the assumption of an angle PHI between the magnetization
vector of surface and tip. The two results are merged in the following way.
For topographies, first create two separate current maps, for ferromagnetic and
antiferromagnetic ordering. Then, move the two CURMAT files to CURFM
and CURAFM, respectively. Now execute bSKAN with the following added
command lines:
MERGE = T
PHI = 45
The angle is given in degree. The result of the calculation is a new file CUR-
MAT containing the merged current maps under this angle. This file can then
be evaluated in the usual manner, extracting current contours, corrugations, or
8.4. MAGNETIC CALCULATIONS 37

current planes. For tunneling spectra the names of the input files to merge are
CURSFM and CURSAFM, apart from that, the procedure is identical.
38 CHAPTER 8. EVALUATION
Chapter 9

FILES

All input and output files are written in capital letters. The first two or three
letters generally specify the information contained. For example, all files related
to the Kohn-Sham states begin with WAV-, files describing the geometry with
A-, input parameters are contained in files named IN-, the output in files OUT-.
Current maps are stored in CUR- files, while plots generally are called PLOT
or PLT.

9.1 Input
Note that we make a difference between essential files (which are necessary for
every run) and non-essential ones (which are usually optional). The input files
are described in Table 13.1 at the end of the guide. The keywords used in
INSCAN file are found in chapter KEYWORDS. For more details on wavefunc-
tions and geometry input files, see chapters WAVEFUNCTIONS and Geometry
files, respectively.

9.2 Output
Output files are either for information on the run, contain the simulation data,
or contain an extracted sample of the simulation data for visualizing. They are
described in Table 13.2

39
40 CHAPTER 9. FILES
Chapter 10

KEYWORDS

The program only uses the first three characters of a keyword, the rest is omitted.
In the following table these essential characters are given in capitals. The list of
essential keywords and their usage can be printed out by executing the program
with only one line in the INSCAN file:
HELP
The current keywords are described in Table 13.3.

10.1 METHOD
10.1.1 TERSOFF-HAMANN
The Tersoff-Hamann method is standard in many DFT simulations, where the
charge within an energy window can be integrated and displayed. The difference
between the bSKAN implementation and DFT implementations is:
1. The surface unit cell is arbitary and can be changed by the keyword CELL
= X Y, which allows in principle to compute an area of multiple unit cells
2. The surface unit cell does not have to be rectangular, while the computed
unit cell is always rectangular. This is an advantage for visualisation
programs like OpenDx
3. The bias dependency can be included from the formulation found for the
first order scattering approach by setting NKELDYSH = 0 (zero order
scattering) or NKELDYSH = 1 (first order scattering).
The line in the INSCAN file has the format:
TERSOFF-HAMANN

10.1.2 STERSOFF
The differential spectrum in this case is calculated with the same model, the bias
dependency for the zero and first order scattering approximation is included via

41
42 CHAPTER 10. KEYWORDS

the NKELDYSH switch. The spectrum can be extended over the whole unit
cell with arbitrary resolution, which allows to compare differential spectra with
locally resolved spectroscopy experiments. The advantages are the same as for
TH spectroscopies.
The line in the INSCAN file has the format:
STERSOFF = V1 V2
here V1 and V2 describe the lower and upper limit of the spectrum. Note
that the bias interval is assumed to contain the zero bias value.

10.1.3 NUMERICAL
Bardeen topographies are based on the wavefunctions of surface and tip; it is
assumed that the z-grid of both systems is equally spaced (0.1 a.u) and that
the distance from the surface nuclei at a given gridpoint i is roughly the same
for surface and tip systems (symmetric setup). The advantage of Bardeen to-
pographies is that they include the electronic structure of the tip explicitly; for
topographies, where typically only a limited number of states around the Fermi
level contributes to the tunneling current, this leads to effects like contrast inver-
sion or contrast changes due to different tip models, a feature well documented
in STM experiments. The bias dependency of the current can be included by
changing NKELDYSH from -1 (the default) to 1 (the first order approximation
in the scattering approach).
The line in the INSCAN file has the format:
NUMERICAL

10.1.4 SPECTROSCOPY
Bardeen spectroscopies include surface and tip electronic structures. The method
is described to some extent in the methods chapter, it is based on differential
increments of the current due to a differential change of the bias voltage. In this
case the dI/dV values are written to a file and then integrated from the point of
zero bias. It is therefore essential that zero bias is included in the calculation.
The binary files CURDSPEC and CURSPEC contain three separate values: one
for the contributions from the surface bandstructure mapped onto the tip Fermi
level, one from the contributions of the tip bandstructure, mapped onto the
Fermi level of the surface, and the sum of the two values.
The line in the INSCAN file has the format:
SPECTROSCOPY = V1 V2
As in the previous cases the bias dependency in the calculation is included
with an appropriate switch NKELDYSH.

10.1.5 FORCE
From version 3.5 the chemical interactions have been included in the simulation
routines. To correct a given CURMAT file for interactions, one needs first to
determine the harmonic constant of surface atoms, and the Wigner Seitz radius
10.2. SETUP 43

of surface and tip atoms. The ratio between current and interaction energy is
parametrized with respect to the Wigner Seitz radii, the parametrization has
the form (see PRL 91, 036803 (2003)):

α = 0.02563 · exp[1.1 · (rS + rT )] (10.1)

The harmonic constants, the distance between nuclei and vacuum boundary of
surface and tip, as well as the lateral position of surface atom and tip apex
are needed to be given in the file INFORCE. At present, the program can only
account for primitive surface cells (one atom only). Then the current values in
CURMAT are used as the basis for a calculation of a file CURMATF, which
contains the currents, corrected for displacement of the surface atoms.
The line in the INSCAN file:
FORCES = T
The default is FORCES = F.

10.1.6 WAVE
Sometimes it is desirable to plot the decay characteristics of a single surface
state. The functionality is provided in the program, in this case the charge
density of a single state, defined by its energy eigenvalue, is plotted for the
on-top position of the unit cell.
The line in the INSCAN file:
WAVE = Ek [htr]
The energy eigenvalue has to be defined precisely enough, so that only a
single state is chosen.

10.2 SETUP
The setup of the calculation involves, apart from the chosen method, the follow-
ing parameters: bias voltage, bias dependency, scan area, thermal broadening,
ferromagnetic and antiferromagnetic transitions, absolute Fermi level, energy
and local resolution of the scan.

10.2.1 ANTIFERROMAGNETIC
In magnetic systems the vector of magnetization will have a direction in space,
which means that spin is no longer isotropic. In this case the calculation of
two separate scans, one with FERROMAGNETIC, the other with ANTIFER-
ROMAGNETIC ordering allows to simulate STM scan on magnetic systems,
where the angle between surface and tip magnetization vectors is used as an in-
put in the subsequent evaluation runs. In non-magnetic systems this parameter
is ignored.
The line in the INSCAN file:
ANTIFERROMAGNETIC
The default in a scan is ferromagnetic ordering of surface and tip states.
44 CHAPTER 10. KEYWORDS

10.2.2 AREA
The tip unit cell is typically made up of a film with a single atomic or a pyramid
apex. This geometry guarantees that the amplitude of the wavefunctions at the
edges of the tip unit cell are negligible compared to the amplitude at the apex.
In this case an integration over the tip unit cell of the overlap of surface and tip
wavefunctions, as required in the Bardeen method, contains mainly the overlap
at the tip apex. The integration area is usually the whole tip unit cell. This is
also the default in every simulation. However, it may be necessary to check the
convergency of an obtained result with respect to the integration area. In this
case the area has to be explicitly specified in units of the tip lattice vectors A1
and A2 . Note that the tip unit cell always has to be at least rectangular.
The line in the INSCAN file:
AREA = a1 a2
The default in every calculation is AREA = 1.0 1.0

10.2.3 BIAS
Setting the bias voltage for a topography simulation, or setting the bias voltage
for evaluation of a spectrum. In the first case the bias voltage either defines the
energy interval for the summation of surface charge (Tersoff-Hamann), or the
shift of Fermi levels of surface and tip systems (Bardeen). In the second case
also the CURRENT value has to be defined; in this case the BIAS/CURRENT
couple defines the setpoint of an STS simulation, as it does in experimental
spectra.
The line in the INSCAN file has the format:
BIAS = Vb
Note that negative bias ranges correspond to tunneling from filled surface
states into empty tip states, positive bias ranges lead to tunneling from empty
tip states into surface states.

10.2.4 CELL
The scan area depends on the symmetry of the surface (see further down) and
the input CELL. For rectangular lattices the variation of the parameter CELL
allows to scan across more than one unit cell. However, since the lateral res-
olution of the scan will be reduced, it is usually more efficient to scan across
a single unit cell and to evaluate the ensuing current matrix over more than
one unit cell. In this case the resolution is retained, while the ensuing constant
current contours still cover a wider area.
The line in the INSCAN file:
CELL = c1 c2
The default in every calculation is CELL = 1.0 1.0
10.2. SETUP 45

10.2.5 DELTA
Defines the broadening σ for the approximation of the delta functional by a
Gaussian in the simulations. The difference between the energy values of surface
and tip states in a scan is calculated and the probability of the transition scaled
with a Gaussian distribution with halfwidth σ. This only applies to topography
simulations, in spectroscopies the energy interval in the differential changes is
commonly small enough (around 20mV) so that this probability is set to one.
The line in the INSCAN file:
DELTA = σ [eV]
The default value in a scan is 100meV.

10.2.6 FERROMAGNETIC
In magnetic systems the vector of magnetization will have a direction in space,
which means that spin is no longer isotropic. In this case the calculation of
two separate scans, one with FERROMAGNETIC, the other with ANTIFER-
ROMAGNETIC ordering allows to simulate STM scan on magnetic systems,
where the angle between surface and tip magnetization vectors is used as an in-
put in the subsequent evaluation runs. In non-magnetic systems this parameter
is ignored.
The line in the INSCAN file:
FERROMAGNETIC
The default in a scan is FERROMAGNETIC ordering of surface and tip
states.

10.2.7 GRIDPOINTS
The lateral resolution in most experimental scans is in the range of 0.1 - 0.2
Å. This means that a typical unit cell of about 3 Å width, can be resolved
by about thirty discrete points of a scan. A simulated scan will scale with the
number of calculated points, so it is generally advisable to limit this number to
the experimentally sensible. The number is defined by the intervals along the
longest axis of the surface scan area. For quadratic unit cells, this leads to the
same number of intervals in both directions, for oblique or hexagonal cells, the
number of intervals in the shorter direction are calculated by the program.
The line in the INSCAN file:
GRIDPOINTS = N
The default is set to GRIDPOINTS = 31. Note that it is possible to do
spectroscopies with only a single gridpoint. In this case the point chosen will
be the point defined by the TOP position (see further down).

10.2.8 HOLLOW
The hollow position in the unit cell in units of the lattice vectors. This cor-
responds usually to the position of surface atoms. The input is used in spec-
46 CHAPTER 10. KEYWORDS

troscopy calculations of a single point. Apart from that it only plays a role in
evaluations.
The line in the INSCAN file:
HOLLOW = t1 t2
The default is set to HOLLOW = 0.5 0.5

10.2.9 LIMITS
Depending on the thermal environment the upper and lower limits of the bias
interval are not sharp but thermally broadened. The keyword LIMITS allows
to adjust the range according to thermal conditions. For room temperature
scans the usual limits will be about 50meV. Note that this parameter is ignored
in spectroscopy, since in this case the differential contributions from one bias
interval to the next is the decisive result of a simulation.
The line in the INSCAN file has the format:
LIMITS = L1 L2
Note also that L1 will be generally negative and L2 positive.

10.2.10 NKELDYSH
The NKELDYSH switch controls whether the scattering approach is used for
calculating the tunneling current, see section on Scattering method. It is highly
recommended to use the zero or first order scattering method for all cases dealing
with non-zero bias, in order to handle the bias dependency correctly.
NKELDYSH = n
where n means the following:

• n = -1 (the default) is the standard Bardeen approach

• n = 0 is the zero order scattering approach

• n = 1 is the first order scattering approach

10.2.11 NSPECTRUM
In spectroscopy runs the bias range is divided in a number of intervals, the dif-
ferential changes are computed separately for every interval. It is recommended
that the energy resolution is about 20-50mV. This is sufficient for tunneling
spectroscopy in the ambient regime, down to about 150 Kelvin. For very low
temperature spectra the surface bandstructure cannot be resolved with sufficient
resolution (required are about 1meV), in this case interpolation routines for the
bandstructure have to be developed. This part of spectroscopy is currently
under development.
The line in the INSCAN file:
NSPECTRUM = N
10.3. EVALUATION 47

The default is set to NSPECTRUM = 11. It should be noted that a lower

resolution does not necessarily increase the speed of the calculation, since tran-
sitions are calculated within every interval. A larger energy interval thus will
increase the number of transitions which have to be calculated.

10.2.12 PIVOT
The keyword is part of a group of three keywords, specifying the scan area. The
PIVOT point is the lower left point of every scan. Its default is PIVOT = -0.5
-0.5, so that the zero point of the surface unit cell is in the middle of the scan.
Note that the units given in PIVOT are units of the lattice vectors A1 and A2
of the surface.
The line in the INSCAN file:
PIVOT = p1 p2
The default in every calculation is PIVOT = -0.5 -0.5

10.2.13 TOP
The on-top position in the unit cell in units of the lattice vectors. This cor-
responds usually to the position of surface atoms. The input is used in spec-
troscopy calculations of a single point. Apart from that it only plays a role in
evaluations.
The line in the INSCAN file:
TOP = t1 t2
The default is set to TOP = 0.0 0.0

10.2.14 ZVACUUM
The vacuum boundary is compulsory to be given for every evaluations in units
of Å as
ZVACUUM = z
The default value is
ZVACUUM = 0.0
which, in turn, results an error message. It should be noted that using
wavefunctions from VASP it is the first value after the STM command in the
INCAR file.

10.3 EVALUATION
bSKAN provides a variety of methods to analyze the data gained in STM/STS
simulations. Generally, it was sought to provide an interface to standard and
open source visualization software. The OpenDX program, which can be down-
loaded free of charge from [Link] is compatible with most of the
output.
The most important command line for an evaluation is
48 CHAPTER 10. KEYWORDS

Figure 10.1: Constant current contour on Si(111) (7 × 7), measured with a

clean tungsten tip. The contour value: 2V, 100pA.

10.3.1 CURRENT
In this case a line can either be:
CURRENT = Ival
or
CURRENT = 0.0
In the first case the evaluation routines search for a current value, which at
the defined BIAS will match Ival . Since the current range in the simulation
depends on the z-range, it is not guaranteed that Ival is part of the result. In
this case the program stops with an errormessage, typically:
CURRENT VALUE NOT FOUND
In this case a line in the OUTSCAN file should specify the minimum and
maximum current value in the simulation. Note that closed current contours
exist only for a limited range of values, so that even though the program writes
a file PLOTCON, the contour may contain holes. The values in the file are
the x, y, and z values of a given Ival contour. They can easily be plotted with
standard utilities, e.g. gnuplot.
In the second case, if Ival = 0.0, the whole three dimensional current map
will be written to a file CURRENT, which contains the current in three nested
loops.
The evaluation routine produces a file which looks exactly like the CHGCAR
files in VASP, with the only differences (i) the z-extension of the lattice, which
is now the z-range of the simulation; (ii) the z-values of the atoms, which are
10.3. EVALUATION 49

now generally negative, since the atoms are below the vacuum boundary. This
file can now either be directly visualized by OpenDX, or the atomic positions
and the current values are combined into three separate OpenDX script files,
with a utility programmed by David Bowler at UCL. In any case the resulting
image looks like (for +2V/50pA on Si(111)) the image shown in Fig. 10.1.

10.3.2 CORRUGATION
It is generally advisable to compute the full current map and to determine the
difference in apparent height for different points on the surface from the plot
routines e.g. OpenDx. However, in simple cases, e.g. on flat metal surfaces,
the corrugation can also be calculated from the difference in apparent height of
two specific points. These points are defined as the TOP and the HOLLOW
position, given in direct lattice coordinates. Here it has to be taken care that
both points are actually part of the surface grid. If they are not, then the nearest
points to the defined ones will be computed by the program automatically. The
input for a corrugation calculation is:
CORRUGATION = T
TOP = xT yT
HOLLOW = xH yH
The output file, PLOTCOR contains apart from the current values at top
and hollow position also the apparent barrier height, defined as the 0.95 the
logarithmic derivative of the current for both positions.

10.3.3 MERGE
For magnetic systems the program calculates the current and spectra depending
on the angle between the surface and tip magnetization vectors. In this case two
separate calculations need to be performed: one with ferromagnetic ordering,
FERROMAGNETIC
and one with antiferromagnetic ordering for the electron transition from
surface to tip. The two current or spectrum maps need to be renamed: the
CURMAT from a ferromagnetic calculation becomes CURFM, the CURMAT
from the antiferromagnetic calculation becomes CURAFM. For spectral maps
the corresponding names are CURSFM and CURSAFM. After the two maps
have been calculated, they can be merged by
MERGE = T
PHI = φ
This creates a new CURMAT or CURSPEC file containing the maps for the
chosen angle φ. This file can then be evaluated in the usual manner

10.3.4 PLOTS
Even though it is generally better to obtain the full current map it is sometimes
necessary to look at the current values at horizontal planes above the surface.
In this case the necessary input is:
50 CHAPTER 10. KEYWORDS

PLOTS = N
FPLOT = FILE
ZPLOT = z1 z2 ... zN
The program then creates N files, the filenames are F ILE.001 to F ILE.N ,
containing the current values in a horizontal plane, specified by the values
ZPLOT.
Chapter 11

Spectroscopy evaluations

A full three dimensional map of all the differential contributions dI(x, y, z, eV )/dV
is probably the most complete information about a surface electronic structure
one can have. From version 3.6 bSKAN is able to compute such a map and to
extract the data in a variety of different manners. From a practical point of
view the information, which can be compared to experimental data is usually:

1. A I(V ) or a dI(v)/dV graph either as a statistical average over a surface

region, or at a specific point of the surface.

2. A two dimensional I(V ) or a dI(V )/dV map over a certain region of the
surface; every gridpoint at the surface given by its unique values.

Experimentally, spectra are always taken at a certain setpoint. This is a

combination of four values: I, V , X, and Y . The I, V values in this case deter-
mine the vertical distance Z. After the setpoint is determined, experimenters
switch off the feedback loop and perform a spectrum I(V ). In certain cases the
bias voltage is oscillated with low frequency and amplitdues of about 20mV.
The dI(V )/dV value is then directly determined by lock-in techniques from the
variation of the current signal.
The three main inputs in an evaluation of a spectrum are the three lines:
BIAS = Vbias
CURRENT = I
TOP = xt yt
The file CURSPEC, which is created after the differential spectrum has been
completed and written to CURDSPEC is then searched for the three input
values. If the combination is not part of the calculation the program stops with
an error message. Otherwise, it performs different tasks, depending on whether
a statistical evaluation is required or not. If the input also contains the line:
STATISTICAL = T
then the whole current and differential current maps are read. The z-value
of the evaluation is then determined depending on the I, V combination in the

51
52 CHAPTER 11. SPECTROSCOPY EVALUATIONS

input. And the output written to PLTSTAT contains the statistical average of
I(V ) and dI(V )/dV values over a two dimensional plane.
The default in evaluations is STATISTICAL = F, in this case the program
writes the separate two dimensional maps to separate files. Here, the convention
is that the dI/dV maps are written to files
[Link]
where xxx is the bias index. The integrated differential currents are written
to files
[Link]
The integrated files are somewhat different than the files one could obtain
by performing a straightforward topography simulation. This is due to the
approximations used in differential spectroscopy. However, the difference should
in general be minor.
In addition, the program writes in both cases a file PLOTSPC, which con-
tains the graph for the spectrum at the TOP position.
Chapter 12

WAVEFUNCTIONS

The format of the wavefunction files used in bSKAN is roughly: (i) a header
giving the scale, the lattice vectors, the number of spins, k-points, maximum
eigenvalues, G-vectors, and z-gridpoints; (ii) for each k-point the position (re-
ciprocal lattice) and the weight, as well as the number of bands and G-vectors;
(iii) the complex amplitudes for each Fourier component at the z-gridpoints in
the vacuum. The first few lines of a typical WAVSAMPLE (WAVTIP) file are
shown in Table 13.4.
The first line is a comment, the next three lines define the surface lattice.
Line one and two give the lattice vectors in the surface plane, line three has only
one element, defining the distance between the first z gridpoint of the amplitudes
and the core of the surface atoms. The next three lines define the electronic
structure of the system: Fermi level, the number of spins, k-points, z-values, G-
vectors (reciprocal lattice vector for the 2dim expansion of the wavefunctions),
and eigenvalues. Then for every k-point the file contains first the coordinates
(in reciprocal space), and the weight of the k-point. Then for every eigenstate
at this point the expansion in reciprocal lattice vectors with their complex and
z dependent amplitudes.

53
54 CHAPTER 12. WAVEFUNCTIONS
Chapter 13

Geometry files

The geometry files ASAMPLE and ATIP contain the atomic positions of the
surface and the tip, respectively in the same format as the CONTCAR file
in VASP simulations ([Link]/vasp/). The first line is a comment
line. The following four lines define the three dimensional repeated unit cell
in the calculation. The next line defines the number of atoms of each species
followed by a line irrelevant for bSKAN, and the final line before the atomic
positions in direct coordinates is:
Direct
It is important to note that cartesian coordinates do not work. It is also
important that the scale (the second line of the file) is equal to unity.
The header of an ASAMPLE (ATIP) file thus looks roughly like:
THIS LINE IS A COMMENT LINE
1.000000000000000
13.5340237999999990 0.0000000000000000 0.0000000000000000
0.0000000000000000 9.5700000000000002 0.0000000000000000
0.0000000000000000 0.0000000000000000 25.0000000000000000
90 2 2
Selective dynamics
Direct
0.1666666666642058 0.0000000000000000 0.0000000000000000 F F F
0.5000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.8333333333357941 0.0000000000000000 0.0000000000000000 F F F

55
56 CHAPTER 13. GEOMETRY FILES
Bibliography

[1] W. Hofer and J. Redinger, Surf. Sci. 447, 51 (2000).

[2] W. A. Hofer, Progs. Surf. Sci. 71, 147 (2003).
[3] W. A. Hofer and A. J. Fisher, Surf. Sci. Lett. 515, L487 (2002).
[4] W. A. Hofer and A. Garcia-Lekue, Phys. Rev. B 71, 085401 (2005).

[5] K. Palotas and W. A. Hofer, J. Phys: Condens. Mat. 17, 2705 (2005).
[6] G. Binnig and H. Rohrer, Helv. Phys. Acta 55, 726 (1982).
[7] G. Binnig, H. Rohrer, C. Gerber, and E. Weibel, Appl. Phys. Lett. 40, 178
(1982).
[8] G. Binnig, H. Rohrer, C. Gerber, and E. Weibel, Phys. Rev. Lett. 49, 57
(1982).
[9] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31, 6207
(1985).
[10] J. Tersoff and D. R. Hamann, Phys. Rev. Lett. 50, 1998 (1985).

[11] J. Tersoff and D. R. Hamann, Phys. Rev. B 31, 805 (1985).

[12] J. C. Chen, Phys. Rev. Lett. 65, 448 (1990).
[13] J. C. Chen, Phys. Rev. B 42, 8841 (1990).
[14] P. Sautet and C. Joachim, Chem. Phys. Lett. 185, 23 (1991).
[15] G. Binnig, N. Garcia, H. Rohrer, J. M. Soler, and F. Flores, Phys. Rev. B
30, 4816 (1984).
[16] F. Tersoff, P. L. de Andres, F. J. Garcia-Vidal, L. Jurczyszyn, N. Mingo,
and R. Perez, Progs. Surf. Sci. 48, 27 (1995).
[17] J. Bardeen, Phys. Rev. Lett. 6, 57 (1961).

57
58 BIBLIOGRAPHY

Table 13.1: INPUT FILES

FILE FUNCTION FORMAT ESSENTIAL
INSCAN all the input parameters ASCII YES
WAVSAMPLE Kohn-Sham states of surface ASCII YES
WAVTIP Kohn-Sham states of model tip ASCII NO
ASAMPLE atomic position of surface atoms ASCII YES
ATIP atomic positions of tip atoms ASCII NO
INEIGENVAL list of eigenvalues excluded ASCII NO
INFORCE list of parameters for FORCE calculation ASCII NO
CURFM current map of ferromagnetic topography ASCII NO
CURAFM current map of antiferromagnetic topography ASCII NO
CURSFM current map of ferromagnetic spectrum ASCII NO
CURSAFM current map of antiferromagnetic spectrum ASCII NO

Table 13.2: OUTPUT FILES

FILE FUNCTION FORMAT ESSENTIAL
OUTSCAN all the output information ASCII YES
CURMAT current matrix of a topographic simulation BINARY YES
CURSPEC current spectrum and matrix of spectroscopy simulation BINARY YES
CURSAVE current spectrum or matrix in ASCII format ASCII YES
TRANSLOG log file for transitions in spectrum simulations ASCII NO
PLOTCON contour plot of topography simulation ASCII NO
PLOTCOR corrugation plot of topography simulation ASCII NO
PLOTSPC I/V spectrum of surface ASCII NO
CURRENT current matrix in OpenDX format ASCII NO
README error messages after input error ASCII NO
BIBLIOGRAPHY 59

Table 13.3: KEYWORDS

KEYWORD FUNCTION
ANTiferromagnetic ordering of surface and tip states
AREa of integration in unit cells of the tip
BIAs voltage in a topography
CELls of the surface in the dimension of the image
CORrugation of the surface electronic structure
CURrent value in the simulated 3D images of the surface
DELta functional mimicked by a Gaussian for tunneling transitions
FERromagnetic ordering of surface and tip states
FORce chemical interactions between surface and tip
FPLot name of the plot files in the output of horizontal plots
GRIdpoints of 2-dim surface structure
HOLlow position on the surface
LIMits in the energy range due to thermal broadening
MERge used for magnetic systems
NKEldysh controls the usage of scattering method
NUMerical evaluation (Bardeen integration)
NSPectrum number of energy gridpoints in a spectrum
PHI angle between magnetization of surface and tip
PIVot point of the surface image
PLOts number of parallel surface plots
STAtistical evaluation of spectra over all surface gridpoints
STErsoff Tersoff-Hamann model of tunneling spectra
SPEctroscopy Bardeen model of tunneling spectra
TERsoff Hamann model of tunneling topographies
TOP position on the surface
WAVe decay characteristics of a single surface state
ZPLot z-gridpoint for the output of parallel plots
ZVAcuum vacuum boundary of the surface in z direction
60 BIBLIOGRAPHY

Table 13.4: WAVFUNCTION FORMAT

ASCII characters
Scale for VASP output: 0.0468651

5.4159547 0.0000000 0.0000000

0.0000000 5.4159547 0.0000000
0.0000000 0.0000000 1.9458969
fermi-energy: 0.0873040
ispin: 2 k-points: 10 z-values: 101 G-vectors: 16 max-eigenval 37
k-point 1 bands 33 G-vectors 13
k-point 0.0625000000 0.0625000000 0.0625000000
eigenenergy 0.0453091754 occupancy 0.0625000000
G-vector: 0 0
( -0.40120E-02, 0.55003E-02)
( -0.40098E-02, 0.54972E-02)
( -0.39832E-02, 0.54609E-02)
( -0.39363E-02, 0.53965E-02)
( -0.38726E-02, 0.53092E-02)
( -0.37954E-02, 0.52033E-02)