Tutorial and Reference Guide: Tomistix OOL IT

ATOMISTIX TOOLKIT
Tutorial and Reference Guide
Version 2.0
Atomistix A/S
Juliane Maries Vej 30
DK–2100 Copenhagen
Denmark
Tel: +45 3532 0630
Fax: +45 3532 0635
Email: [email protected]
Atomistixr and Virtual NanoLabr are registered trademarks of Atomistix A/S.

Atomistix ToolKitTM is a trademark of Atomistix A/S.
Copyright c 2003–2006 Atomistix A/S.

All rights reserved.
Revision date: April 24, 2006 (2.0.3)

Preface
The strength of the Atomistix ToolKit (ATK) software lies in its flexibility to describe sys-
tems with different symmetries. It can describe isolated systems (molecules), periodic
systems (crystals), and systems of the type bulk–nanodevice–bulk (two-probe systems).
In this tutorial, we will first start with a quick tour, which outlines the capabilities of
the software, and then focus on each of the different symmetries, with emphasis on the
two-probe geometry.
The ATK software uses a novel design which eases the setup of two-probe calculations.
However, these types of calculations add complexity to the electronic structure calcu-
lations, and a range of input parameters are required. The purpose of this document is
to give the user a hands-on introduction to the ATK package, and to provide a guide
to the various commands and input parameters. A detailed summary of all input and
output parameters and available commands can be found in the appendices.
The background material, describing the underlying physics and algorithms, is de-
scribed in a number of research papers listed in the References section (page 125)
of the manual. The electron transport part is mainly described in Refs. 1, 11, 13. In this
tutorial we will not delve deep into the physics details, but focus on describing how to
use the program.
The proper citation to ATK includes a reference to the most important methodology
papers. We recommend users to cite:
ATK version 2.0, Atomistix A/S (www.atomistix.com).
M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, Phys. Rev. B 65,
165401 (2002)
J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and

D. Sánchez-Portal, J. Phys. Condens. Matter 14, 2745 (2002).
J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 245407 (2001).
Remember to visit our website http://www.atomistix.com for updates to this manual

and the software!
Atomistix, April 2006.

Contents
1 Introduction 1
1.1 About This Tutorial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Quick Tour 5
2.1 Molecular Orbitals of an H2 Molecule . . . . . . . . . . . . . . . . . . . 6
2.2 Band Structure of a Lithium Chain . . . . . . . . . . . . . . . . . . . . . 7
2.3 Setting up the Two-Probe System . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Electrode Calculation . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.2 Determining the Li–H Bond Length . . . . . . . . . . . . . . . . 12
2.3.3 Placing H2 Between the Lithium Chains . . . . . . . . . . . . . . 13
2.3.4 Self-Consistent Calculation at Zero Bias . . . . . . . . . . . . . . 14
2.4 Analyzing the Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5.1 MPSH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.2 Density of States . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.6 Current–Voltage Characteristics . . . . . . . . . . . . . . . . . . . . . . 17
2.6.1 Voltage Drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.7 Relaxation of a Two-Probe System . . . . . . . . . . . . . . . . . . . . . 19
2.8 Including a Gate Electrode . . . . . . . . . . . . . . . . . . . . . . . . . 20
3 Molecular Systems 23
3.1 Water Molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 The Self-Consistent Loop . . . . . . . . . . . . . . . . . . . . . . . . . . 24

iv Contents
3.2.1 How To Cure Convergence Problems . . . . . . . . . . . . . . . 25
3.3 Relaxations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 Changing the Accuracy of the Calculation . . . . . . . . . . . . . . . . . 27
3.5 Molecular Unit Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.6 Spin-Polarized Calculations . . . . . . . . . . . . . . . . . . . . . . . . 30
4 Periodic Systems 35
4.1 Band Structure of a Si Crystal . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Convergence of an Al Wire . . . . . . . . . . . . . . . . . . . . . . . . . 38
5 Two-Probe Systems 41
5.1 Geometry of an Ideal Lithium BCC Lattice . . . . . . . . . . . . . . . . 41
5.1.1 Lithium Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.1.2 Two-Probe Geometry . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.3 Miswrapped Electrode Atoms . . . . . . . . . . . . . . . . . . . 45
5.1.4 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . 45
5.1.5 Energy Integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.2 Benzene Ring Coupled to Lithium Electrodes . . . . . . . . . . . . . . . 46
5.2.1 Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2.2 Central Region . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.2.3 Finite Bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3 Heterogeneous Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . 51
6 Spin-Polarization and k-point Sampling of Two-Probe Systems 59
6.1 Transmission Through a Fe–MgO–Fe Interface . . . . . . . . . . . . . . 59
6.1.1 Fe Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.1.2 Fe-MgO-Fe Geometry . . . . . . . . . . . . . . . . . . . . . . . 61
6.1.3 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . 62
6.2 k-point Sampling in Two-Probe Systems . . . . . . . . . . . . . . . . . . 62
6.3 k-point Sampling of the Transmission Spectrum . . . . . . . . . . . . . . 65
6.3.1 Plotting the k-dependent Transmission . . . . . . . . . . . . . . 65

Contents v
Appendix A ATK vs. TranSIESTA-C 69
A.1 Changes to Keywords . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
A.1.1 Keywords No Longer Available . . . . . . . . . . . . . . . . . . 71
A.1.2 Keywords With New Names . . . . . . . . . . . . . . . . . . . . 72
A.1.3 Keywords With New Names and Syntax . . . . . . . . . . . . . 72
Appendix B Basis Sets 75
B.1 Basis Set Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
B.1.1 Single Zeta Orbitals . . . . . . . . . . . . . . . . . . . . . . . . . 75
B.1.2 Double Zeta Orbitals . . . . . . . . . . . . . . . . . . . . . . . . 76
B.1.3 Polarization Orbitals . . . . . . . . . . . . . . . . . . . . . . . . 76
B.2 Pseudopotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
Appendix C Command Line Options 79
C.1 Input and Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
C.2 Environment Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Appendix D Parameters and Keywords 83
D.1 Syntax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
D.1.1 Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
D.1.2 Units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
D.2 Keywords . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
D.2.1 Analysis:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
D.2.2 AtomList:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
D.2.3 Bulk:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
D.2.4 Electrode:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
D.2.5 Method:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
D.2.6 Molecule:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
D.2.7 NumOrb:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
D.2.8 Relaxation:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
D.2.9 SCF:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
vi Contents
D.2.10 Simulation:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
D.2.11 System:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
D.2.12 TwoProbe:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
D.2.13 UnitCell:: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Appendix E Output file 105
E.1 Header and self-consistent cycle . . . . . . . . . . . . . . . . . . . . . . 105
E.2 Total energy and atom list . . . . . . . . . . . . . . . . . . . . . . . . . 107
E.3 Mulliken population . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
E.4 Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
E.4.1 Energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . 111
E.4.2 Transmission spectrum . . . . . . . . . . . . . . . . . . . . . . . 111
Appendix F Writing Scripts 113
F.1 Bash Scripts for Linux and Cygwin . . . . . . . . . . . . . . . . . . . . . 113
F.2 Batch Files for Windows . . . . . . . . . . . . . . . . . . . . . . . . . . 115
F.3 Python Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Appendix G Atomic Data 121
References 125
Index 127
C H A P T E R
1
INTRODUCTION
ATOMISTIX TOOLKIT is a library of atomic scale modeling techniques that can be used to
calculate a wide range of properties of nanoscale systems. New methods are constantly
being added to the toolkit, and for each new release ATK will be able to calculate more
properties using a wider range of different modeling techniques. The most unique
features is the ability to calculate the electrical properties of nanoscale devices, which
consist of a scattering region coupled to two macroscopic bulk systems or electrodes.
Atomistix provides two different interfaces for performing electronic transport calcula-
tions: VIRTUAL NANOLAB (VNL), which is a graphical user interface (GUI), and ATOM-
ISTIX TOOL K IT (ATK), described in this manual, which is a file-based or command line
interface. Which of the interfaces to use, is sometimes a matter of taste, and sometimes
determined by the task at hand. Most users will probably find that VNL is the most ap-
pealing option, because of its intuitive ease of use. Some prefer ATK because it is more
transparent, gives more control, and calculations can be automated by using scripts.
Moreover, ATK provides immediate access to the newest functionality, whereas VNL
instruments for accessing these new features are often developed in a later step.
A modeling study generally consists of three parts: setting up the calculation, running
the calculation and finally analyzing the output of the calculation. The purpose of this
tutorial is to help the user getting control of these steps using ATK, often in combination
with VNL.
When setting up a calculation, the first thing to consider is the symmetry of the system.
ATK provides methods for describing three different system types; molecular, bulk, and
two-probe systems. A two-probe system consists of a nanoscale region coupled to two
macroscopic bulk systems, or electrodes. The nanoscale region can be a molecule,
nanotube, cluster of atoms, a piece of a semiconductor or an interface between the
electrodes. In this tutorial, we will first discuss the modeling of molecular and periodic
1
2 Introduction
Figure 1.1: ATK is capable of modeling isolated, periodic, and open systems using both very
accurate first principle modeling methods and very fast semi-empirical methods.
systems before introducing the more complex two-probe systems.
Once the system to be studied is defined, the calculational methodology must be spec-
ified. The description of the electronic structure in ATK is based on density functional
theory (DFT) [5] and an important ingredient in the calculations is the approximation
used for the exchange–correlation functionals.
In the current version of ATK, the local density approximation (LDA) [8] and the gener-
alized gradient approximation (GGA), in the form of the Perdew–Burke–Ernzerhof (PBE)
exchange–correlation functional [7], are implemented. The analogue spin-dependent
exchange–correlation functional are implemented as well, which make it possible to
perform spin-polarized calculations. For most applications, the PBE functional will be
the more accurate option.
Other sources of errors arise from the use of pseudopotentials to describe the effect
of the core electrons, the finite size of the basis set and the numerical integration
routines. These errors can in most cases be systematically reduced at the expense
of the computational time. The default parameters for ATK are chosen such that the
approximation for the exchange–correlation potential usually is the largest source of
error. In this tutorial, we will describe the other sources of errors, and explain how
the default parameters can be changed in order to reduce the computational cost or
increase the accuracy.
1.1 About This Tutorial 3
Figure 1.2: A typical two-probe system, consisting of a carbon nanotube between two metal
surfaces. The unique capability of ATK is its ability to study the electron transport through such
a system.
Finally, the output of the calculation must be analyzed to gain new insight from the
streams of data produced by the program. ATK can be used to calculate many different
properties of the system, including electronic current, voltage drop, transmission coef-
ficients, electron density, etc. The results can then be exported to standard file formats
so that they can be visualized in two or three dimensions, e.g. in VIRTUAL NANOLAB.
1.1 ABOUT THIS TUTORIAL
The first part of this tutorial is a quick tour, outlining the capabilities of ATK. In the
following three sections, we will go into more details about the available commands
and parameters relevant for different system symmetries. If you are new to ATK, it is
recommended to read the tutorial in its entirety, since many basic concepts will be
introduced and explained in the quick tour and then used in the subsequent sections.
We shall assume that the software has been properly installed and the environment
variables set (including the path) as described in the Installation Guide, supplied with
the program. The Installation Guide is also available from Atomistix’ web site.
We recommend that you follow the tutorial, type the input file and enter the commands
manually to gain experience with the steps needed for an ATK calculation. It is also
recommended to create a new directory where you will run the examples. However, for
reference, all input files and scripts used in this tutorial are located in the subdirectories
of the directory examples in the program installation directory.
In this document, we symbolize the command prompt by a dollar sign $. Thus, state-
ments like
$ atk h2.atk h2.out
are commands to be entered at the command prompt. In Linux and similar envi-
ronments a command can be broken into multiple lines by using the backslash (’\’)
character, cf. the box on page 114. In Windows this is not possible, and it is necessary
to type the entire command line on a single line, remembering to omit the backslash.
4 Introduction
For more detailed information about command line and input file parameters, please
consult the reference sections, Appendices C and D.
Note that commands in the input file are not case sensitive. The capitalizations used
throughout this manual are only meant to ease the reading.
Note! When editing files, be sure to save them as plain text files,
using ANSI or ASCII encoding.
C H A P T E R
2
QUICK TOUR
In this first part of the tutorial, you will be guided through a few examples that illustrate
the general features and capabilities of ATK. The focus will be on physical properties
that you can study with ATK, rather than the particulars of all commands used to cal-
culate them. Nevertheless, we shall try to explain the role and influence of most of the
parameters used in the calculations when they appear.
The aim of the calculations is to study the conductance of a two-probe system con-
sisting of a hydrogen molecule coupled with two lithium chains, as illustrated in Fig-
ure 2.1. Our motivation for considering this structure is a recent measurement of the
conductance of a single hydrogen molecule [10]. In the experiment, the molecule was
situated between two platinum wires, but for understanding the system, it is possible
to use lithium chains instead, which reduces the computational requirements substan-
tially. We will investigate, why this system has a nearly perfect conductance despite
the huge hydrogen HOMO–LUMO gap.
To get experience it is recommended to type in the input file for these calculations your-
self, but they can also be found in the directory examples/quicktour in the program
installation directory.
Figure 2.1: A hydrogen molecule between two semi-infinite one-dimensional lithium chains.
5
6 Quick Tour
2.1 MOLECULAR ORBITALS OF AN H2 MOLECULE
Create a text file h2.atk with the following content
System::Type Molecule
Simulation::Type Relaxation
AtomList::Format Angstrom
%block Molecule::AtomList
H 0.0 0.0 -0.35
H 0.0 0.0 0.35
%endblock Molecule::AtomList
and execute the command
$ atk h2.atk h2.out
Note that no output is displayed on the screen while the program is running (except
for two lines indicating that the program started and finished), since we have specified
a second file, h2.out, to which all output is redirected. If this second parameter is
omitted, the corresponding output is instead displayed on the screen. Hence, you may
alternatively redirect or pipe the output to another file, system device, or program by
use of the standard operators ’>’ and ’|’.
The first line of the input file tells ATK that we wish to study a molecular system and the
second line instructs the program to calculate the equilibrium positions of the hydrogen
atoms (i.e. to relax the geometry). This is achieved through a number of self-consistent
total energy and force calculations where the position of the hydrogen atoms are moved
along the force directions. The initial guess for the distance between the atoms is set
to 0.7 Å, as seen from the Molecule::AtomList block. The coordinates are given in
Ångström (10–10 m), as indicated by the keyword AtomList::Format.
Note! If you mistyped any of the parameters in the input file,

the program will most likely fail to run, and will produce an er-
ror message specifying the offending input. Correct the typos and
try again. Note that by default ATK will not overwrite old output
files, and will refuse to run unless you either delete the old output
files or specifically tell it to overwrite them by using the --force
option.
Two output files are generated from the calculation:
h2.out Contains a log of the calculation and important results such as total energies
can be extracted from this file.
2.2 Band Structure of a Lithium Chain 7
h2.nc NETCDF data file with the final state of the calculation and data for plotting
with e.g. Atomistix’ visualization program VIRTUAL NANOLAB (see below). The
file contains a complete description of the final state of the calculation and can
thus be used to restart a calculation or to generate further analysis data. For more
information about NETCDF files, see Appendix C.1.
We can now find the relaxed distance between the two hydrogen atoms by looking for
the word "Positions" in the output file h2.out. The obtained result is approximately
0.77 Å, which is in acceptable agreement with the experimental value of 0.741 Å, in
particular in view of the fact that the hydrogen molecule is the smallest molecule of all,
and therefore poorly described by the local density approximation (LDA), which is the
default exchange–correlation functional in ATK.
ATK can also be used to extract further results from the calculation by specifying ad-
ditional analysis options. Keywords starting with Analysis:: do not affect the actual
electronic structure calculation, but only produce analysis output that is stored in a
NETCDF file. This information can then be visualized in VIRTUAL NANOLAB.
If we, for example, wish to inspect the orbitals of the HOMO and LUMO levels, we
create an input file h2_analysis.atk
Analysis::MolecularOrbital::A::Index 1
Analysis::MolecularOrbital::A::Label Lumo
Analysis::MolecularOrbital::B::Index 0
Analysis::MolecularOrbital::B::Label Homo
and run ATK with the command
$ atk -a h2.nc h2_analysis.atk h2_analysis.out
The option -a tells the program to read the data file h2.nc where the results of the elec-
tronic structure calculation are stored, and use the data from this file for the analysis.
The molecular orbitals are written to the file h2_analysis.nc, and this file can now be
imported into VIRTUAL NANOLAB to plot the eigenfunctions. Such a plot is illustrated in
Figure 2.2. Note that the electron density and electrostatic potential are always written
to the NETCDF file by default.
2.2 BAND STRUCTURE OF A LITHIUM CHAIN
In the previous section, we looked at the first ingredient we need for our final two-
probe system, namely the hydrogen molecule that will form the central region. Now
we shall consider the leads, which in our example are semi-infinite one-dimensional
chains of lithium atoms.
8 Quick Tour
Figure 2.2: An isosurface of the LUMO orbital of the hydrogen molecule, created using Atomistix
visualization program VIRTUAL NANOLAB.
The following input file (lichain.atk) defines a lithium chain.
System::Type Bulk
NumOrb::MeshCutoff 100 Ry
NumOrb::BasisSet::Size SZ
Bulk::NumKPoints::C 200
UnitCell::LatticeConstant 3.0 Angstrom
%block Bulk::UnitCell
3. 0. 0.
0. 3. 0.
0. 0. 1.
%endblock Bulk::UnitCell
%block Bulk::AtomList
Li 0.0 0.0 0.0
%endblock Bulk::AtomList
Analysis::BandLine::L::Start (0., 0., 0.0) FracRecip
Analysis::BandLine::L::End (0., 0., 0.5) FracRecip
Analysis::BandLine::L::NumPoints 1000
In this case, we have a bulk system as indicated by the first line in the input file. The
single Li atom is repeated periodically in three dimensions according to the lattice vec-
tors defined by the unit cell. In ATK, the three directions defined by the lattice vectors
are referred to as the A, B and C directions, respectively. This notation is reflected
in a number of keywords, such as Bulk::NumKPoints::C used above to specify the
number of points that should be used in the direction of the reciprocal lattice vector
corresponding to the C direction.
In the input file above, the lattice vectors defining the unit cell are given in units of
the lattice constant, defined by UnitCell::LatticeConstant and thus, the lithium
2.2 Band Structure of a Lithium Chain 9
−10.102
−10.103
Total Energy (eV)
−10.104
−10.105
−10.106
2.82 2.84 2.86 2.88 2.9 2.92
Lattice Constant (Å)
Figure 2.3: Total energy versus the lattice constant of the lithium chain.
atom is repeated by 9, 9, and 3 Å in the x, y, and z directions, respectively. Lithium

atoms separated by 9 Å are too distant to interact significantly (the lattice constant 3.0 Å
roughly corresponds to the experimental value) and hence, we have defined an infinite
chain of lithium atoms extending in the z direction (or rather the C direction, but these
coincide in the present case).
The command
$ atk lichain.atk lichain.out
calculates the total energy of the lithium chain for a lattice constant of 3.0 Å.
Inspecting the output file, one notes that it contains substantially more parameters than
we specified in the input file. Input parameters which are not specified explicitly in the
input file are set to default values specified in the file $ATK_DATA_DIR/atk-defaults.
atk. The default values are also listed in Appendix D.2.
For example, if the parameters NumOrb::MeshCutoff and NumOrb::BasisSet::Size
are removed from lichain.atk, ATK will use the default parameters. This will actually
result in more accurate results; however, the calculation will also take more time to
execute. In particular, the SZ basis set is too small to yield an even reasonably accurate
description of the band structure of lithium. It is nevertheless often useful to start with
a small basis set to get a general picture of the system, before launching a larger and
more accurate, but also more time-consuming calculation.
In the next step, we will use ATK to find the equilibrium lattice constant of the chain,
i.e. the value of the lattice constant that minimizes the total energy. In this case, it is
not possible to use the relaxation method employed for the hydrogen molecule, since
this approach only relaxes the atoms within the unit cell and not the size of the unit
10 Quick Tour
5
2.84 Å
2.90 Å
4
3
Energy (eV)
−1
−2
0 0.2 0.4 0.6 0.8 1 1.2
kz (1/Å)
Figure 2.4: Band structure of a lithium chain with an interatomic distance of 2.84 Å and 2.90 Å,
respectively. The Fermi level is at zero energy, and thus the Li chain is metallic.
cell itself.
Instead, we employ the fact that ATK takes a number of command line options, as
we have already seen examples of. Of particular interest for our task is the option
--param "parameter". Using this when calling ATK is equivalent to including the
same statement "parameter" in the input file. In this way, it becomes possible to write
small scripts that perform complex tasks.
The scripts can be created in any suitable scripting language, but here we will fo-
cus on the ATK commands, rather than the implementation of the scripts. However,
in Appendix F, a collection of example scripts are provided for Linux/Cygwin (bash
shell scripts), Windows (DOS batch files) and platform-independent scripts written in
Python, for each of the tasks discussed in the tutorial.
To calculate the total energy for a given lattice constant (2.84 Å in this case), we use
the command
$ atk lichain.atk lichain2.84.out -o lichain2.84.nc \

--param "UnitCell::LatticeConstant 2.84 Angstrom"
The command may be written on a single line (without the ’\’), it is just broken here to
be easier to read. We are already familiar with the first line, which contains the names
of the input and output files. The option -o is required to create NETCDF files with
unique names for each lattice constant (the default name is the same as the input file),
and finally we specify the value of the lattice constant using the --param option, as
discussed above.
A script can be written (see Appendix F) that runs such a calculation for a sequence of
2.3 Setting up the Two-Probe System 11
lattice constants in the range from 2.82 to 2.92 Å. In this way, a plot like the one in
Figure 2.3 can be made. The figure shows that the chain with lattice constant 2.90 Å has
the lowest energy and thus, we will use this lattice constant for the lithium electrodes in
the two-probe system. Note, that the total energy in itself is meaningless, only relative
total energies are useful.
The last three lines in the script lichain.atk computes the band structure of the
lithium chain between the point, defined by k=(0, 0, 0), and the X point, defined
by k = (0, 0, /a), where a is the lattice constant (or rather, in our case, the spacing of
the atoms in the C direction). The keyword FracRecip means that the wave vectors,
such as (0.0, 0.0, 0.5), are interpreted in units of the reciprocal lattice vectors. Other
units are also possible, cf. Appendix D.1.2.
The band structure is shown in Figure 2.4 for two different lattice constants; note that
the calculations only need to be carried out for positive wave numbers kz , since the
band structure always is symmetric about the point.
2.3 SETTING UP THE TWO-PROBE SYSTEM
Our goal is to place the hydrogen molecule between two lithium chains and calculate
the transmission through the system. Thus, the system consists of three parts, the left
electrode, the central region, and the right electrode. The left and right electrodes are
semi-infinite lithium chains, while the central region consists of the hydrogen molecule
with parts of the lithium chains. The geometry is illustrated in Figure 2.5.
2.3.1 ELECTRODE CALCULATION
We will first consider the electrode atoms, i.e. the red atoms in Figure 2.5. For this we
need to do a calculation with three-dimensional boundary conditions of the electrode
atoms, similar to the lithium chain calculation in the previous section. However, for
an electrode calculation, it is assumed that there are only interactions between atoms
in neighboring unit cells, meaning that all other interactions are set identically to zero.
If, on the other hand, the physical interactions in reality have a longer reach than
Figure 2.5: Geometry of the hydrogen molecule coupled with semi-infinite lithium chains. The
electrode atoms are colored red.
12 Quick Tour
the lattice constant (this is the situation that would occur in most realistic cases), it
is necessary to create a conventional cell, or supercell, such that all interactions are
properly accounted for.
In the present case, we may safely assume that lithium atoms separated by more than 4
lattice constants are not interacting significantly. Therefore, we include 4 lithium atoms
in the supercell (see the input file below), which then of course has a four times larger
lattice constant in the C direction. In this way, 4 atoms in one unit cell may interact
with the 4 atoms in an adjacent cell, and we thus account for all interactions extending
up to at least 4 of the (original) lattice constants along the chain.
To define the electrodes, create a file lielec.atk with the content
System::Type Electrode
Electrode::NumKPoints::C 500
%block Electrode::UnitCell
3. 0. 0.
0. 3. 0.
0. 0. 4.
%endblock Electrode::UnitCell
%block Electrode::AtomList
Li 4.35 4.35 0.00
Li 4.35 4.35 2.90
Li 4.35 4.35 5.80
Li 4.35 4.35 8.70
%endblock Electrode::AtomList
and produce the electrode NETCDF file lielec.nc with the command
$ atk lielec.atk lielec.out
2.3.2 DETERMINING THE LI–H BOND LENGTH
The next step is to place the hydrogen molecule between the lithium chains, and we
therefore need to estimate the Li–H bond length. We do this by relaxing a Li4 H2 cluster.
For this purpose, create a file li4h2.atk with the content
Relaxation::ForceTolerance 5.e-2 eV/Ang
Relaxation::MaxSteps 200
2.3 Setting up the Two-Probe System 13
NumOrb::BasisSet::Size::Li SZ
Li 0.00 0.00 -4.955
Li 0.00 0.00 -1.905
H 0.00 0.00 -0.375
H 0.00 0.00 0.375
Li 0.00 0.00 1.905
Li 0.00 0.00 4.955
Here we introduced some new parameters related to relaxation. They determine the
tolerance (in the forces) for deciding when the calculations are considered to have
converged, and how many steps this at most may take. We also specified a different
(smaller) basis set specifically for the lithium atoms by appending the element symbol,
Li, to the keyword NumOrb::BasisSet::Size. In this way, we still use the larger default
basis set for the hydrogen atoms, which results in a more accurate calculation of the
relaxed H–H bond length.
To perform the calculation, execute the command
$ atk li4h2.atk li4h2.out
The calculation takes a while to complete.
From the output file, we find that the relaxed Li–H distance is estimated to be 2.37 Å,
and the H–H bond is slightly stretched to 0.80 Å. In the next section, we use the values
to set up the central region of the two-probe system.
2.3.3 PLACING H2 BETWEEN THE LITHIUM CHAINS
We now place the hydrogen molecule between the lithium contacts. In the central
region, we have H2 coupled with three lithium atoms to each side. For the Li–H–H–Li,
part we will use the relaxed geometry from the Li4 H2 cluster, which can be extracted
from the output file li4h2.out. The remaining Li–Li distances are set to the lithium
chain lattice constant found earlier (2.90 Å). The file lih2-0.0.atk defines the input
to ATK for the complete two-probe system:
System::Type TwoProbe
Method::Type NumOrb
NumOrb::BasisSet::Size::Li SZ
14 Quick Tour
%block TwoProbe::CentralAtomList
Li 0.0 0.0 -8.568
Li 0.0 0.0 -5.668
Li 0.0 0.0 -2.768
H 0.0 0.0 -0.402
H 0.0 0.0 0.402
Li 0.0 0.0 2.768
Li 0.0 0.0 5.668
Li 0.0 0.0 8.568
%endblock TwoProbe::CentralAtomList
TwoProbe::LeftElectrode::NetCDFFile lielec.nc
TwoProbe::RightElectrode::NetCDFFile lielec.nc
TwoProbe::LeftSurface::NumberOfAtoms 3
TwoProbe::RightSurface::NumberOfAtoms 3
The unit cell is automatically generated by ATK from the electrode file lielec.nc and
there is therefore no need to specify the unit cell in the input file. The electrodes are
automatically aligned with the central region using the electrode lattice constant.
To check that the system is set up correctly we will generate a data file with the two-
probe geometry. With the -c option ATK does not perform an actual calculation. Type
the command
$ atk -c lih2-0.0.atk --xyz
The --xyz option is used to output the geometry into lih2-0.0.xyz, which can be
visualized with the Atomistix program VIRTUAL NANOLAB (see Figure 2.5).
2.3.4 SELF-CONSISTENT CALCULATION AT ZERO BIAS
The structure looks correct and next we perform the self-consistent calculation:
$ atk lih2-0.0.atk lih2-0.0.out --force
Note that this calculation will take rather long time. Note that we have included the op-
tion --force to instruct ATK to overwrite and replace any existing output and NETCDF
files. This is necessary in this case, since the previous command (at the end of the
previous section) already produced an output NETCDF file with the same name as the
one we will create here.
It is instructive to inspect the output file. The first part of the file contains information
about the system setup and input parameters. Below this, there is information about
the self-consistent cycle. Two self-consistent cycles are performed. The first cycle is
2.4 Analyzing the Results 15
carried out with periodic boundary conditions, and the only purpose of this part is to
obtain a good initial guess for the density matrix. The second cycle is carried out with
the proper open system boundary conditions.
2.4 ANALYZING THE RESULTS
As already mentioned above, ATK offers a range of analysis options, all specified
by Analysis:: keywords, which can be used to post-process the results of the self-
consistent energy structure calculation. The data generated can be used to draw con-
clusions about the results, and also for plotting.
We shall analyze the transmission spectrum of the system, and for this we create a file
trans.atk with the contents
Analysis::TransmissionSpectrum::NumPoints 100
Analysis::TransmissionSpectrum::E0 -5. eV
Analysis::TransmissionSpectrum::E1 5. eV
Analysis::MPSH::A::Index 0
Analysis::MPSH::A::Label Mpsh0
Analysis::MPSH::B::Index 1
Analysis::MPSH::B::Label Mpsh1
Analysis::DOS::A::Energy 0 eV
Analysis::DOS::A::Label total
Analysis::DOS::A::Type total
and run the command
$ atk trans.atk -a lih2-0.0.nc trans.out
The results will be saved in trans.nc, and some results can also be inspected in the
log file trans.out. Note how the output file name for the NETCDF file is created from
the input atk file, unless a specific name is specified by using the option -o.
2.5 TRANSMISSION SPECTRUM
The main quantity of interest in most two-probe calculations is the transmission spec-
trum. This is controlled by the keywords in trans.atk beginning with Analysis::
TransmissionSpectrum. The resulting spectrum is plotted in Figure 2.6. We observe,
as indicated earlier, that the transmission is very high (above 0.6) at the Fermi level.
This is a great puzzle, since the HOMO–LUMO gap of a hydrogen molecule is huge
compared to typical electron excitation energies.
16 Quick Tour
1.0
0.8
Transmission Coefficient
0.6
0.4
0.2
0.0
−5 −4 −3 −2 −1 0 1 2 3 4 5
Energy (eV)
Figure 2.6: Transmission coefficients of an H2 molecule placed between two lithium chains. The
dashed line shows the Fermi level.
2.5.1 MPSH
We shall now try to use the other analysis options specified in trans.atk to explain
the observed unusually high conductance of the hydrogen molecule in the lithium
chain. A very useful concept is the molecular projected self-consistent Hamiltonian
(MPSH). The molecular projection refers to the part of the two-probe system which
remains when the electrode and surface atoms are excluded. In the input file lih2-0.
0.atk, we defined three surface atoms on each side of the hydrogen atoms by using the
options TwoProbe::LeftSurface::NumberOfAtoms and TwoProbe::RightSurface::
NumberOfAtoms. These two parameters specify the number of surface atoms in the
TwoProbe::CentralAtomList. If their values are n and m, respectively, the first n
atoms in the central atom list will be treated as surface atoms of the left electrode and
the last m atoms in the list will be treated as surface atoms of the right electrode.
Thus in our case, all lithium atoms are either electrode or surface atoms and the molec-
ular projection only refers to the hydrogen molecule. This means that the MPSH is the
Figure 2.7: The two small red and turquoise ellipsoids in the center of the image show an
isosurface (the color denotes the phase) of the wave function of the MPSH state with energy
2.36 eV. The phase difference clearly indicates that this state corresponds to the LUMO. The
large balls are not part of the isosurface, but illustrate the positions of the lithium atoms.
2.6 Current–Voltage Characteristics 17
self-consistent Hamiltonian of the hydrogen molecule in the presence of the lithium

chain. The new eigenvalues are –8.02 eV for the HOMO and 2.36 eV for the LUMO.
The reduction of the HOMO–LUMO gap is directly responsible for the large value of
the transmission coefficient of the hydrogen molecule observed experimentally [10].
In the file trans.atk, we included the lines Analysis::MPSH::A::Index and Analysis::

MPSH::A::Label. These lines mean that the MPSH wavefunctions are written in the
analysis file trans.nc and can then be visualized with VIRTUAL NANOLAB. In Fig-
ure 2.7, we show a plot of the MPSH state with energy 2.36 eV, the MPSH LUMO
orbital.
2.5.2 DENSITY OF STATES
From the previous section, we conclude that the transmission is mainly due to the
LUMO, since the position of this orbital is closest to the transmission maximum at
around 2.5 eV (cf. see Figure 2.6). Another piece of evidence to support the conclusion
that the transmission mainly takes place through the LUMO is obtained from the density
of states (DOS) at the Fermi level. The last three lines in trans.atk, starting with
Analysis::DOS, control whether or not the DOS is output to the file trans.nc. The
corresponding data has been visualized, spatially resolved, with VIRTUAL NANOLAB in
Figure 2.8 , and we see that there is a substantial density of states concentrated to the
hydrogen atoms at the Fermi level.
2.6 CURRENT–VOLTAGE CHARACTERISTICS
The next task is to calculate the current–voltage characteristics of the system. This re-
quires a self-consistent calculation for a series of different values of the applied bias,
and for this purpose we may again create a script (Appendix F), with the central com-
mand (assuming a bias of 0.8 V)
$ atk lih2-0.0.atk lih2-0.8.out -o lih2-0.8.nc -i lih2-0.6.nc \

--param "TwoProbe::LeftElectrode::Voltage 0.8 eV"
Figure 2.8: Volume plot of the density of states of the LiH2 system at the Fermi energy (0 eV).
The small spheres indicate the positions of the atoms.
18 Quick Tour
40
Single Configuration
Relaxation
30
Current (µA)
20
10
0
0 0.2 0.4 0.6 0.8 1
Bias (V)
Figure 2.9: I–V curve of a hydrogen molecule between two lithium chains, comparing the
current–voltage characteristics when the atomic positions are relaxed (cf. Section 2.7) or not.
Note carefully the -i option, which refers to a NETCDF file for a different bias (0.6 V
in this case), which we thus assume has already been calculated (cf. the scripts). In this
way, i.e. by using the -i option, we can employ a previously calculated density matrix
as the initial guess for the new calculation to improve the calculation speed.
The command above adds an additional parameter to the lih2-0.0.atk parameter

file. The parameter TwoProbe::LeftElectrode::Voltage adds a bias voltage to the
left electrode. The resulting current can then be extracted from the output file. Look
for the line containing the word "Current". A positive value of the current indicates that
the current flows from the left electrode to the right; obviously the electron flow is in
the opposite direction.
Figure 2.10: Voltage drop through the Li–H–H–Li system for an applied bias of 1 V.
2.7 Relaxation of a Two-Probe System 19
0.820 2.38
H−H
H−Li
Li−H
H−Li, Li−H Bond Length (Å)

H−H Bond Length (Å)
0.815 2.36
0.810 2.34
0 0.2 0.4 0.6 0.8 1
Bias (V)
Figure 2.11: Change in the H–H and Li–H bond lengths as a function of the bias.
2.6.1 VOLTAGE DROP
Next, we calculate the voltage drop through the structure. The voltage drop is given by
the difference in the effective potential between the finite bias calculation and the zero
bias calculation. The effective potential of the zero bias calculation was stored in the
file trans.nc.
We extract the effective potential from the 1.0 V calculation with the command
$ atk trans.atk trans+1.0.out -o trans+1.0.nc -a lih2-1.0.nc
In Figure 2.10, the difference between the effective potentials from the two calculations
are plotted. Note, that for this effectively one-dimensional system the central region
(the number of surface lithium atoms) is not large enough to screen the applied electric
field, and therefore an abrupt change occurs in the voltage drop at the right electrode.
2.7 RELAXATION OF A TWO-PROBE SYSTEM
In this section, we relax the two-probe geometry both at zero and finite bias.
Make a new file lih2-0.0-relax.atk with the same content as lih2-0.0.atk, and
add the line
20 Quick Tour
5 HOMO (Self−Consistent)
LUMO (Self−Consistent)
HOMO (Non Self−Consistent)
LUMO (Non Self−Consistent)
0
MPSH Eigenvalues (eV)
−5
−10
−15
−20
−10 −8 −6 −4 −2 0
Gate Voltage (eV)
Figure 2.12: The curved lines represent the self-consistent change in the MPSH eigenstates as
function of applied gate voltage, while the straight lines are linear extrapolations (i.e. the initial
shift of the eigenstates by the gate voltage) of the MPSH eigenstates at zero gate voltage.
We can reuse the script we used to calculate the current–voltage relationship, and just
replace the input file with the one containing the relaxation command above. Thus,
the generic ATK command becomes (for the bias 0.8 V, using the 0.6 V as initial guess)
$ atk lih2-0.0-relax.atk lih2-0.8-relax.out \

--force -o lih2-0.8-relax.nc -i lih2-0.6-relax.nc\
--param "TwoProbe::LeftElectrode::Voltage 0.8 eV"
These calculations take a few hours on a modern PC. If you check the output file,
you will find that only the hydrogen atoms have moved (the results are summarized
in Figure 2.11). This is because the positions of the electrode and surface atoms are
kept fixed. Remember that surface atoms are specified by the commands TwoProbe::
LeftSurface::NumberOfAtoms and TwoProbe::RightSurface::NumberOfAtoms.
2.8 INCLUDING A GATE ELECTRODE
As the final step of the quick tour, we demonstrate how to simulate the electrostatic
potential induced by a third electrode (a gate electrode) on the two-probe system.
The strength of the electrostatic potential due to the gate is controlled by the keyword
TwoProbe::Molecule::GateVoltage.
First, create a copy of the input file lih2-0.0.nc named gate-0.0.nc. This file will
be used as the initial guess for the lowest gate voltage. We can now write a script
2.8 Including a Gate Electrode 21
(see Appendix F), which performs a self-consistent calculation of the electronic levels
for our two-probe system for gate voltages ranging from –1 V to –10 V by using the
command (here at the gate voltage –5 V)
$ atk lih2-0.0.atk gate-5.0.out -o gate-5.0.nc -i gate-4.0.nc \

--param "TwoProbe::Molecule::GateVoltage -5.0 eV"
The gate voltage for ATK should be specified as an electrostatic potential energy (here
in electron volts). The conversion of the gate voltage used in a physical system to this
potential energy is not trivial. It requires a specification of the geometry used in the
physical system and a solution of the Poisson equation that relates the electrostatic
potential on the molecule to the applied gate voltage.
The gate electrode is not included as a physical electrode in ATK, and, consequently,
a current cannot run from the source or drain electrode to the gate electrode. In the
current implementation, we simulate the electrostatic effect of the gate electrode by
simply shifting the MPSH part of the Hamiltonian with the gate voltage (i.e. the electro-
static potential energy). This corresponds to assuming that the gate electrode induces
an external potential localized to the molecular region. For metallic electrodes this will
usually be a reasonable approximation.
Figure 2.12 shows the self-consistent shift of the MPSH eigenstates due to the applied
gate potential. For reference is also plotted the initial (non-self-consistent) shift. We see
that the self-consistent change in the MPSH eigenvalues is significantly smaller than
the initial shift by the gate potential.
22 Quick Tour
C H A P T E R
3
MOLECULAR SYSTEMS
In this chapter, we go into more details about how to use ATK to study molecular sys-
tems. The purpose of the different examples is to show the syntax of the input file and
to demonstrate how different input parameters affect the accuracy of the calculations.
We recommend that you type the input files yourself, but they can also be found in the
directory examples/molecule in the program installation directory.
3.1 WATER MOLECULE
Make an input file h2o-1.atk for an H2 O molecule with the following lines
Method::Type NumOrb
H 0.757 0.586 0.000
H -0.757 0.586 0.000
O 0.000 0.000 0.000
Since the input file does not explicitly list parameters related to the basis set or which
method to use for the calculation, the default values will be used. These correspond
to using density functional theory (DFT) with the LDA exchange–correlation functional
of Perdew and Zunger [8], and a double zeta basis with polarization orbitals for all the
elements. Run the file
23
24 Molecular Systems
$ atk h2o-1.atk h2o-1.out
and inspect the results in the output file.
3.2 THE SELF-CONSISTENT LOOP
An important part of the DFT calculation is to obtain the self-consistent density, and a
number of parameters are dedicated to controlling this part of the calculation. We will
now take a look at some of the parameters. Add the following two lines to the previous
input file
SCF::Coordinate DensityMatrix
SCF::DiagonalMixingParameter 0.5
and run ATK on the new input file, which we call h2o-2.atk
By inspecting the output file h2o-2.out, you will find that ATK now converged in fewer
steps.
The parameter SCF::Coordinate controls which system variable is mixed during the
self-consistent cycle. For closed systems with a fixed number of electrons, i.e. molec-
ular or bulk geometries, usually density matrix mixing is most efficient. For open two-
probe systems, where the central region does not have a fixed number of electrons,
Hamiltonian mixing is, based on our experience, most efficient.
The algorithm used for mixing the system variable is controlled by the keyword (SCF::
Algorithm). In ATK, there are two different mixing algorithms, Broyden [2] and Pulay
mixing [9]. Both schemes use a number of previous iterations to make a guess on the
new density. The mixing ratio between the new guess and the old guess is controlled by
SCF::DiagonalMixingParameter. If the value of this parameter is too large, the self-
consistent loop will not converge. If it is too small the self-consistent loop will take a
very long time. For isolating systems, large values of SCF::DiagonalMixingParameter
can be used (up to 0.5), whereas for metallic systems the value should usually be
around 0.1.
All parameters which control the self-consistent cycle are summarized in Table 3.1.
3.3 Relaxations 25
SCF::Coordinate DensityMatrix k Hamiltonian k EffectivePotential

This keyword determines which coordinate is used for mixing.
SCF::Algorithm Pulay k Broyden
This keyword controls the mixing algorithm.
SCF::Tolerance real
When the relative change of the density and the Hamiltonian matrix is less than this
value between two iterations, the self-consistent cycle is considered converged.
SCF::DiagonalMixingParameter real
This parameter defines the fraction of the new value of the SCF::Coordinate which
is mixed with the old value.
SCF::HistorySteps (integer)
Number of previous steps used in the mixing algorithm to generate the next state.
SCF::MaxSteps integer
Maximum number of steps in the self-consistency loop.
Table 3.1: Parameters used to control the self-consistent cycle. For units and default values, refer
to the keyword reference section on pp. 86– 104.
3.2.1 HOW TO CURE CONVERGENCE PROBLEMS
Convergence problems with molecular systems can usually be related to a small or

vanishing HOMO–LUMO gap (the HOMO is highest occupied molecular orbital, the
LUMO (the lowest unoccupied molecular orbital). In such cases, a very small perturba-
tion can shift a molecular state from occupied to unoccupied, making the convergence
difficult. To remove this instability one can increase the electron temperature. Usually
a value of
Molecule::ElectronTemperature 0.2 eV
will be sufficient to obtain convergence.
If an increased electron temperature does not cure an convergence problem, increasing

the number of history steps used in the Pulay mixing to e.g. 20 (default is 6), might
work. The number of history steps is specified by the keyword SCF::HistorySteps.
Moreover, lowering the value of SCF::DiagonalMixingParameter can also improve
the convergence properties of a system; the required value may be as low as 0.01.
3.3 RELAXATIONS
ATK calculates the forces on all atoms, and this information can be used to relax the
ionic coordinates, as we already demonstrated in the quick tour. The relaxation module
can be used for molecules, bulk or two-probe geometries.
Create a file h2o-3.atk with the parameters of h2o-1.atk, and add the following lines:
Relaxation::ForceTolerance 5.e-2 eV/Ang
Relaxation::MaxSteps 200
Now the command
will produce a relaxed geometry.
It is possible to fix the position of some of the atoms during the relaxation using the key-
word Simulation::ConstrainAtoms. The following command will fix the positions
of the two hydrogen atoms (atoms number 0 and 1) during the relaxation:
$ atk h2o-3.atk h2o-3c.out -o h20-3c.nc \

--param "Simulation::ConstrainAtoms 0:1"
Relaxation::ForceTolerance real ForceUnit

The relaxation is completed when the largest force component is less than this value. The
sX j
total force magnitude is also returned in the output file. This is the length of the force
vector on all atoms, i.e.
fi;k j2
i;k
where fi;k is the force component k on atom i, and i = 0; :::; n 1 where n is the number
of atoms and k = x; y; z .
Relaxation::MaxSteps integer
The maximum number of energy evaluations during the relaxation run.
Relaxation::MaxDisplacement real LengthUnit
Sets the maximal displacement for each step in the relaxation. If this is too large the
relaxation might fail, while a too small step will slow the relaxation down.
Simulation::ConstrainAtoms int:int[,int:int]
Fixes the positions of certain atoms during a relation. Example: a system consists of 12
atoms numbered 0, 1, 2, : : :, 11 (the numbering of atoms always begins at zero). With the
command
Simulation::ConstrainAtoms 0:2,7:10,4:5
we constrain atoms 0, 1, 2, 7, 8, 9, 10, 4, 5. Negative atom indices correspond to counting
backwards from zero, so the command
Simulation::ConstrainAtoms 0:2,-5:-2,-8,5
is equivalent to the one above.
Table 3.2: Parameters used to control the relaxation. For units and default values, refer to the
keyword reference section on pp. 86–104.
3.4 Changing the Accuracy of the Calculation 27
3.4 CHANGING THE ACCURACY OF THE CALCULATION
Usually the main approximation in an ATK calculation is related to the choice of the
exchange–correlation functional. The exchange-correlation potential is chosen using
the keyword NumOrb::XCFunctional. Currently four different exchange–correlation
functionals are available.
LDA-PZ The local density approximation (LDA) with the Perdew–Zunger parametriza-
tion [8] of the correlation energy of a non spin-polarized homogeneous electron
gas calculated by Ceperly–Alder [4].
LSDA-PZ The local spin density approximation (LSDA) with the Perdew–Zunger parametriza-
tion [8] of the correlation energy of a spin-polarized homogeneous electron gas
calculated by Ceperly–Alder [4].
GGA-PBE The non spin-polarized generalized gradient approximation (GGA) in the

Perdew–Burke–Ernzerhof approximation (PBE) [7].
SGGA-PBE The spin-polarized generalized gradient approximation (SGGA) in the Perdew–

Burke–Ernzerhof (PBE) approximation [7].
In the previous section, we relaxed the H2 O molecule with the LDA-PZ functional.
Inspection of the relaxed geometry shows that the H2 O molecule has a bonding an-
gle of HOH =103.3 and an OH bond length of d OH =0.97 Å, which are close to the
experimental values of HOH =104.5 and d OH =0.96 Å. We will now test the GGA-PBE
functional (the spin-polarized functionals, LSDA-PZ and SGGA-PBE will be described in
Section 3.6). To change the exchange–correlation functional add the following line
NumOrb::XCFunctional GGA-PBE
to the h2o-3.atk and save it in a new file h2o-4.atk
The new values HOH =103.0 and d OH =0.97 Å are not much different from the previous
ones. In general, structural properties do not depend much on the exchange-correlation
functional (i.e. LDA versus GGA, or LSDA versus SGGA). In contrast, energy changes
are more sensitive to the choice of exchange-correlation functionals, and the GGA-PBE
functional will usually perform better than the LDA functional.
The accuracy of the calculation is also affected by the choice of the pseudopotential,
the basis set and the numerical integration routine, and it is important to check that the
errors due to these parameters are small. In this section, we will show how the size of
the basis set and the accuracy of the numerical integration is controlled.
To change the basis set to a minimal basis set, add NumOrb::BasisSet::Size SZ to

the input file h2o-3.atk and save the file as h2o-5.atk. Relax the H2 O molecule
defined by this file. By inspection of the output file (h2o-5.out) we find that the
geometry has changed. For instance, HOH =99.1 , whereas it was HOH =103.3 for the
DZP basis set. Thus, the use of the SZ basis set decreases the accuracy of the calculation
substantially, whereas DZP basis sets usually are comparable to well-converged plane-
wave basis sets.
ATK has 5 different built-in basis sets:
SZ A single basis orbital for each valence orbital.
DZ Two basis orbitals for each valence orbital.
SZP SZ and one basis orbital for the first unoccupied shell.
DZP DZ and one basis orbital for the first unoccupied shell.
DZDP DZ and two basis orbitals for the first unoccupied shell.
It is possible to use different basis sets for different atoms. For instance, the following
keywords (h2o-6.atk)
NumOrb::BasisSet::Size DZ
NumOrb::BasisSet::Size::O DZP
assign a DZP basis to the oxygen atoms and a DZ basis to all other atoms. It is possible
to control the basis sets in much more detail, in order to improve the results further;
this is described in Appendix B.
The parameter NumOrb::MeshCutoff controls the density of the grid used for real space
integrals. The real space resolution, x is related to the mesh cutoff, Ecut through the
p
relation x = = Ecut , where Ecut is in Rydberg and x in bohr.
Now make the following calculation
$ atk h2o-1.atk h2o-1-50Ry.out -o h2o-1-50Ry.nc \

--param "NumOrb::MeshCutoff 50 Ry"
You will note that the calculation is much faster now, but the total energy has increased.
The cutoff value of 50 Ry is far too low to give any meaningful results for the total
energy. The default value of NumOrb::MeshCutoff is 150 Ry. Let us try to increase the
value to 200 Ry:
$ atk h2o-1.atk h2o-1-200Ry.out -o h2o-1-200Ry.nc \

--param "NumOrb::MeshCutoff 200 Ry"
In Figure 3.1, we show the total energy and the force in the water molecule as a
function of the mesh cutoff. The force converges more slowly than the total energy. For
most materials the system default of 150 Ry will be sufficient.
3.5 Molecular Unit Cell 29
−465.25 1.05
Largest force component (eV/Å)

−465.30
Total energy (eV)
1.00
−465.35
Total Energy
Force
−465.40 0.95
0 50 100 150 200 250 300
NumOrb::MeshCutoff (Ry)
Figure 3.1: The total energy of H2 O and the largest force component (which always is on the
oxygen atom, in the z -direction) as a function of the mesh cutoff. The geometry has not been
relaxed.
3.5 MOLECULAR UNIT CELL
In ATK, a molecule is embedded in a unit cell with periodic boundary conditions. This
means that even molecular systems are treated as periodic systems, and there is in fact
not only a single molecule in the system, but an infinite number of molecules.
The size of the unit cell, and thus the distance between the repeated copies of the
molecule, is controlled by the parameter Molecule::ElectrostaticPaddingFactor.
When this parameter is set to zero, the unit cell has the minimal size where there are
no matrix elements, i.e. overlap between wave functions of molecules in neighboring
cells, which in turn means that they only interact electrostatically. The default value is
0.1.
Try to increase the size of the unit cell with the command
$ atk h2o-1.atk h2o-1-pad0.7.out -o h2o-1-pad0.7.nc \

--param "Molecule::ElectrostaticPaddingFactor 0.7"
By inspecting the output file you will find that the total potential energy has only
changed slightly, whereas the molecular eigenvalues are shifted by about 0.04 eV.
In Figure 3.2, we show the dependence of the total energy and the HOMO energy
as a function of the padding factor. The change in the computed values is due to an
electrostatic interaction between the repeated images of the molecules.
We see that the total energy is well converged with the default choice of padding factor.
The H2 O molecule has a relatively strong dipole moment, so the default parameter
−465.32 −6.03
−6.04
HOMO Energy (eV)

−6.05
Total Energy (eV)
−465.33 −6.06
−6.07
−6.08
Total Energy
HOMO
−465.34 −6.09
0 0.2 0.4 0.6 0.8 1 1.2
Molecule::ElectrostaticPaddingFactor
Figure 3.2: The total energy of H2 O as a function of the value of Molecule: :

ElectroStaticPaddingFactor , and the energy of the HOMO. Note that the scale on the two
y-axes is the same.
should be sufficient for most systems. The HOMO energy shows a larger variation with
the padding factor. This is because the molecular eigenvalues are reported relative to
the vacuum level (the value of the effective potential far away from the molecule), and
the determination of the vacuum level is more sensitive to the electrostatic interaction
between the cells.
Thus, for molecules with a large dipole moment the padding factor must be increased
to obtain converged molecular eigenvalues relative to the vacuum level.
3.6 SPIN-POLARIZED CALCULATIONS
In this section, we introduce calculations with spin dependent exchange-correlation

functionals by considering the O2 molecule.
First, make the file o2.atk with the context
AtomList::Format Ang
O 0.0 0.0 -0.6
O 0.0 0.0 0.6
NumOrb::XCFunctional SGGA-PBE
3.6 Spin-Polarized Calculations 31
$ atk o2.atk o2.out
By inspecting the output file, we see that there are two sets of eigenvalues, one for spin
up and one for spin down, and that there are two more electrons with spin up than spin
down. Hence, the electronic ground state is a triplet in agreement with experimental
evidence (or, at least, the spin along the quantization axis is MS = 2 12 ~).
In the spin dependent calculations, the initial density is generated by superposition

of atomic densities with maximum spin. The calculation then automatically relaxes
the spin to that of the ground state, although for complicated systems, it might be
trapped in a local minimum. We can modify the behaviour by using the keywords
SCF::InitialSpinPolarization and Molecule::TotalSpin, cf. Appendix D.2 The
former one specifies the initial spin, whereas the latter constrains the total spin to a
given value. Below, we constrain the total spin in order to study the O2 molecule in
the quintet state by copying o2.atk to o2_5.atk and adding the line
Molecule::TotalSpin 2 hbar
Note, that the total spin is given in units of ~ and that the spin of a single electron is 12 ~.
Hence, the above line indicates that there should be 4 electrons in excess with spin up,
the majority spin, corresponding to a quintet state. Next, we execute the command
$ atk o2_5.atk o2_5.out
An inspection of the output file reveals that the molecule is stretched significantly com-
pared to the previous calculation without a constraint on the total spin.
The bond lengths from the two calculations and the corresponding energies are listed
in Table 3.3 together with the results from a calculation of the singlet state, obtained
by executing
$ atk o2.atk o2_1.out -o o2_1.nc --param "Molecule::TotalSpin 0 hbar"
and a calculation with the spin-independent GGA exchange-correlation functional

(GGA-PBE), obtained by executing
$ atk o2.atk o2-ns.out -o o2-ns.nc --param "NumOrb::XCFunctional GGA-PBE"
As can be seen in Table 3.3, the triplet state is indeed the groundstate, whereas the
quintet and singlet states are not stable. Moreover, computing the singlet state with a
spin-dependent exchange-correlation functional really just corresponds to performing
a spin-indepedent calculation as already indicated in Table 3.3. However, the spin-
independent calculation can be made at lower cost.
State XC Distance (Å) Total Energy (eV)

Singlet GGA-PBE 1.22 –871.50
Singlet SGGA-PBE 1.22 –871.50
Triplet SGGA-PBE 1.22 –872.71
Quintet SGGA-PBE 1.91 –867.25
Table 3.3: Relaxed O2 molecule in different spin states and computed with spin-dependent or
independent exchange-correlation functionals.
For some of the Analysis keywords, it is now necessary to specify the spin of interest.
For instance, if we want to compare the HOMO of the triplet structure with that of the
quintet structure, we may write an input file (o2-homo3.atk)
Analysis::MolecularOrbital::A::Label HOMO
Analysis::MolecularOrbital::A::Spin Up
$ atk o2-homo3.atk o2-homo3.out -a o2.nc
in order to get the HOMO of the triplet structure. To get the HOMO of the quintet
structure, we instead use (o2-homo5.atk)
Analysis::MolecularOrbital::A::Label HOMO
Analysis::MolecularOrbital::A::Spin Up
$ atk o2-homo5.atk o2-homo5.out -a o2_5.nc
Surface plots of the two orbitals are plotted in Figure 3.3. It is seen that the HOMO
of the triplet state is an anti-bonding orbital, whereas that of the quintet state is an
Figure 3.3: The HOMO of the triplet (left panel) and quintet (right panel) states of O2 illustrated
with VIRTUAL NANOLAB. Note that the O–O distance is 1.2 Å in the left panel, whereas it is
1.9 Å in the right panel.
3.6 Spin-Polarized Calculations 33
anti-bonding orbital. In fact, by studying the eigenvalues and plotting the corre-
sponding orbitals, it can be realized that in the quintet state a spin-down electron in a
degenerated and bonding orbital (the spin down eigenvalues with index 3 and 4) has
been transferred to an anti-bonding orbital (the HOMO in the quintet state). This in
turn explains the increased inter-atomic distance in the quintet state.
In general, spin-polarized calculations come at a higher cost, but in many cases they
are necessary in order to get meaningful energies and structures as indicated by the
above example. Other examples are individual atoms, radicals, and some systems con-
taining transition metals like iron. An example of the latter is the spin filter described
on pp. 59–67.
C H A P T E R
4
PERIODIC SYSTEMS
In a periodic system, the eigenfunctions obey Bloch’s theorem and can be labeled by
the wavevector k. In this section, we will show how k-point sampling is implemented
in ATK, and how the number of k-points can be increased to obtain well converged
results.
4.1 BAND STRUCTURE OF A SI CRYSTAL
Type the following lines into the file si-1.atk (the input files for this chapter can be
found in examples/bulk in the installation directory)
System::Type Bulk
Numorb::XCFunctional GGA-PBE
Numorb::basisset::Size SZ
# FCC LATTICE
UnitCell::LatticeConstant 5.43 Ang
0. 0.5 0.5
0.5 0. 0.5
0.5 0.5 0.
Bulk::NumKPoints::A 3
Bulk::NumKPoints::B 3
35
36 Periodic Systems
# TWO POINT BASIS

AtomList::Format Scaled
AtomList::Scale 5.43 Angstrom
Si 0.00 0.00 0.00
Si 0.25 0.25 0.25
This input file demonstrates an optional way of specifying the atomic coordinates in
the atom list. By using the command AtomList::Format Scaled, we express all coor-
dinates in units of the length AtomList::Scale, which here is set to 5.43 Å. Also note
that the unit cell is defined in terms of the lattice constant.
Also note the lines starting with the symbol #. This is the comment symbol, and the
line will not be read by the input parser.
Run the file with the command
$ atk si-1.atk si-1.out
The command performs a self-consistent calculation for the silicon crystal with 3 3
3=27 k-points uniformly spaced in the first Brillouin zone (i.e. in the reciprocal lattice)
following the Monkhorst–Pack scheme [6]. The three directions defined by the lattice
vectors, specified in the unit cell block in the input, file are known in ATK as the A, B
and C directions, respectively. The present case is a good example showing that these
directions do not always coincide with the Cartesian axes (which of course are used to
define the components of the lattice vectors). The three lines containing the keyword
NumKPoints specify the density of the k-point sampling in the directions of the three
reciprocal lattice vectors corresponding to the A, B, and C directions, respectively.
If a, b, and c are the lattice vectors, then the reciprocal lattice vectors ka , kb and kc are
defined as
ka = 2(b c)=V; kb = 2(c a)=V; kc = 2(a b)=V;
where V = |a (b c)| is the volume of the unit cell. In the case of Si, which has a
diamond crystal structure corresponding to fcc (face-centered cubic) with two atoms in
the basis, the reciprocal unit cell is the same as the real space unit cell of bcc (body-
centered cubic), and thus
ka = 2a ( 1; 1; 1); kb = 2a (1; 1; 1); kc = 2a (1; 1; 1);
where a is the Si lattice constant as given in the input file above.
We would now like to plot the band structure from the point k = (0, 0, 0) to the X point
k = (2 /a, 0, 0). Clearly, the X point corresponds to k = (kb + kc )/2, a relation which
is independent of the actual value of the lattice constant. Therefore it is convenient to
use the unit FracRecip to specify this point.
4.1 Band Structure of a Si Crystal 37
10
0
Energy (eV)
−5
−10
−15
SZ−GGA
DZP−GGA
−20
Γ X
k (Å)
Figure 4.1: The band structure of a Si crystal along the –X direction calculated with GGA.
To calculate the band structure along the line –X, create the file siband.atk with the
contents
System::Type Bulk
UnitCell::LatticeConstant 5.43 Ang

0. 0.5 0.5
0.5 0. 0.5
0.5 0.5 0.
Analysis::BandLine::A::Start (0., 0., 0.) FracRecip

Analysis::BandLine::A::End (0., 0.5, 0.5) FracRecip
Analysis::BandLine::A::NumPoints 100
and run the command
$ atk siband.atk si-1-band.dat -a si-1.nc -o si-1-band.nc
It is crucial to inform ATK about the unit cell when using the FracRecip unit, otherwise
the reciprocal lattice vectors cannot be calculated correctly. If the unit cell is omitted
from the input file in analysis mode, the program will use a default unit cell corre-
sponding to a simple cubic system with lattice constant 1 Å. However, the program
will not read the unit cell from the input file unless we also state that the system type is
bulk (this is done by the first line in the input file). If one does not use the FracRecip
unit, but instead specifies the start and end points of the bandline explicitly in inverse
Bohr or Angstrom, then one may exclude the unit cell from the input file.
38 Periodic Systems
Figure 4.1 shows the obtained Si band structure, which indicates that Si would have an
indirect band gap of 1.73 eV. Experimentally, however, Si is known to have an indirect
band gap of 1.13 eV. The reason for the discrepancy is, in this case, that the SZ basis
set is not versatile enough to describe a Si crystal well enough. The SZ basis set is
spherically symmetric, but a Si crystal is held together by strongly directed covalent
bonds. Therefore, the next thing to try is to use a DZP basis set, by the command:
$ atk si-1.atk si-2.out -o si-2.nc --param "NumOrb::BasisSet::Size DZP"
Then, extract the band structure by the command
$ atk siband.atk si-2-band.dat -a si-2.nc -o si-2-band.nc
With the more accurate basis set we find a band gap below 0.6 eV, i.e. the band gap
is severely underestimated. This is a well known artifact of DFT, or at least LDA and
GGA-PBE.
4.2 CONVERGENCE OF AN AL WIRE
In this section, we will calculate the total energy of a one-dimensional aluminium wire
and illustrate how the k-point sampling can be increased to converge the total energy.
Create the following file (alchain-1.atk)
System::Type Bulk
Method::Type NumOrb
Bulk::ElectronTemperature 1.e-3 Ry
UnitCell::LatticeConstant 1. Ang
15. 0. 0.
0. 15. 0.
0. 0. 2.8
Al 7.5 7.5 1.4
and run the calculation
$ atk alchain-1.atk alchain-1.out

4.2 Convergence of an Al Wire 39
−140.05
−140.10
Total Energy (eV)
−140.15
−140.20
−140.25 8
12
20
40
−140.30 60
120
0.01 0.1
Bulk::ElectronTemperature (eV)
Figure 4.2: Total energy as function of the electronic temperature. For each electron temperature
there are data for runs with 8, 12, 20, 40, 60, and 120 k-points.
We will now check if the total energy is converged with respect to the number of k-
points in the chain direction (the C direction). Try to perform runs with 8, 12, 20, 40,
60 and 120 k-points by changing the value of the keyword Bulk::NumKpoints::C.
Use the command
$ atk alchain-1.atk alchain-2.out -o alchain-2.nc \

--param "Bulk::NumKPoints::C 20"
or write a script to perform the whole series of calculations (cf. Appendix F).
In Figure 4.2, we show the total energy as a function of the number of k-points and
the electron temperature, and we see that more than 20 k-points are needed to ob-
tain a well converged result at the low electron temperature used in the calculation
(alchain-1.atk, 10–3 Ry = 0.013 eV).
The problem with the Al chain is that it is a metal and the Fermi level is therefore
crossing a band. The region around the Fermi level needs to be described accurately,
which requires a dense k-point sampling. By introducing a finite electron temperature,
the Fermi level crossing is smoothed out and fewer k-points are needed.. The elec-
tron temperature is changed with the command Bulk::ElectronTemperature. The
following command generates a run with 10 k-points and an electron temperature of
0.1 eV.
$ atk alchain-1.atk alchain-t0.1.out -o alchain-t0.1.nc \

--param "Bulk::ElectronTemperature 0.1 eV"
ATK accepts energies given as temperature-equivalents, and it is often convenient to

specify the electron temperature as e.g. 300 K.
40 Periodic Systems
Relative Energy (eV) 0.10
0.01
8
12
20
40
60
120
0.00
0.01 0.10
Bulk::ElectronTemperature (eV)
Figure 4.3: Total energy per atom of an Al chain where the bond lengths are alternatingly 2.7 Å
and 2.9 Å, relative (on an absolute scale) to the energy of an Al chain where all bond lengths are
equal to 2.8 Å (cf. Figure 4.2). Results are shown for calculations using 8, 12, 20, 40, 60 and
120 k-points (symbols as in Figure 4.2). Note the different energy scales in Figure 4.2 and in this
figure.
In Figure 4.2, the total energy is plotted as function of the electron temperature for
different numbers of k-points. As can be seen, the total energy is a function of both
these quantities, and for a given temperature it is important to choose sufficiently many
k-points to get a converged result, otherwise, artifacts may arise. However, as indicated
by Figure 4.2, a higher electron temperature facilitates convergence of the total energy
with respect to the k-points.
Relative energies converge faster with respect to number of k-points than total energies.
Moreover, increasing the electron temperature will typically change a relative energy
only slightly, which in turn can be used to reduce the number of k-points further. This
is illustrated in Figure 4.3, where it can be seen that the relative energy is less sensitive
to the choice of k-points than the total energy (compare the different energy scales of
Figures 4.2 and 4.3), and that the relative energy only depends slightly on the electron
temperature.
C H A P T E R
5
TWO-PROBE SYSTEMS
The unique capability of ATK is the calculation of so-called two-probe systems, and we
shall in this last part of the tutorial devote special attention to setting up and analyzing
these nanoscale systems.
A two-probe geometry consists of two semi-infinite systems coupled via a scattering

region. In ATK, this geometry is broken down into three subsystems: the left electrode,
the central region, and the right electrode. Figure 5.1 illustrates this division of the
geometry for a lithium bcc crystal in the (100) direction. The central region must
include the part of the Li (100) surfaces perturbed by the scattering region (indicated
in the figure by a dashed line). Beyond the dashed line, the perturbation from the
scattering region is screened out by the metallic leads, and the effective potential retains
its bulk value. Thus, the effective potential of the complete system can be obtained from
two independent calculations of the electrode regions, and a composite two-probe
calculation of the central region.
5.1 GEOMETRY OF AN IDEAL LITHIUM BCC LATTICE
In the following, we describe the setup of the system illustrated in Figure 5.1. The input
files are, for reference, located in examples/tutorial in the installation directory.
41
42 Two-Probe Systems
5.1.1 LITHIUM ELECTRODES
The Li(100) electrodes are defined by the following input file (li-100.atk)
Method::Type NumOrb
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 4.0
AtomList::Format Fractional
Li 0.0 0.0 0.0
Li 0.5 0.5 0.125
Li 0.0 0.0 0.25
Li 0.5 0.5 0.375
Li 0.0 0.0 0.5
Li 0.5 0.5 0.625
Li 0.0 0.0 0.75
Li 0.5 0.5 0.875
We here introduce yet another way of specifying the atomic coordinates in the atom
list; the Fractional keyword implies that all coordinates are expressed in terms of the
lattice vectors, defined in the UnitCell block. As before, we use the simplest basis set
(SZ) in order to speed up the calculations in these examples. For more accurate results,
Figure 5.1: Two-probe geometry of a Li (100) crystal. The left electrode, the
central region, and the right electrode are combined by matching the hatched
atoms with the electrode lattice vectors. The hatched atoms are defined
through the keywords TwoProbe::LeftElectrode::EquivalentElectrodeAtom and
TwoProbe::LeftElectrode::EquivalentCentralAtom.
5.1 Geometry of an Ideal Lithium BCC Lattice 43
use the default value by excluding the line which specifies the basis set from the input
file.
Lithium is a bcc metal, and in the (100) direction there are two inequivalent planes, A
and B. The electrode file sets up the layer structure (ABABABAB). There is one atom per
layer, with periodic boundary conditions.
In choosing the electrode cell, three important considerations must be taken.
The cell must be long enough in the transport direction, such that matrix elements
only extend to the first repeated cell in the z direction. This is because ATK
assumes that all other matrix elements are zero; if they, in reality, are not, the
contribution from these overlaps will not be properly accounted for.
The cell must be periodic in all three directions.
The C direction must be perpendicular to the A and B directions.
These requirements, in particular that of a periodic electrode cell in the C direction in

conjunction with the first point, may sometimes make it necessary to use a very long
electrode cell. However, this only increases the computational efforts of generating
the electrode, while the two-probe calculation is not affected, since the two-probe
calculation only includes orbitals in the electrode region that are interacting with the
central region.
We next calculate the effective potential of the electrode region with the command:
$ atk li-100.atk li-100.out
Since the electrode region of the left and the right electrode is the same, we can use
the same electrode file for both.
5.1.2 TWO-PROBE GEOMETRY
We next set up the central region with the commands (li-ideal.atk)
Method::Type NumOrb
# ----------------------------------
# Ideal Li (100) system
# ----------------------------------
TwoProbe::LeftElectrode::NetCDFFile li-100.nc
TwoProbe::RightElectrode::NetCDFFile li-100.nc
TwoProbe::Contour::DeltaEMin 2.00000 Ry
AtomList::Format Ang
Li 0.00000 0.00000 0.00000
Li 1.75500 1.75500 1.75500
Li 0.00000 0.00000 3.51000
Li 1.75500 1.75500 5.26500
Li 0.00000 0.00000 7.02000
Li 1.75500 1.75500 8.77500
Note that no unit cell parameters are included in the two-probe setup file; this infor-
mation, including the k-point sampling, is obtained from the electrode files.
The next step is to align the central region with the electrode atoms. ATK makes
an automatic alignment as illustrated in Figure 5.1. The automatic alignment re-
lies upon the specification of equivalent atoms in the electrode and in the central
region. These atoms are defined with the keywords TwoProbe::LeftElectrode::
EquivalentCentralAtom and TwoProbe::LeftElectrode::EquivalentElectrodeAtom.
The two equivalent atoms will be positioned one electrode lattice vector apart.
By default the first atom in the central region is positioned one lattice vector from the
first atom in the left electrode cell, and the last atom in the central region is positioned
one lattice vector from the last atom in the right electrode cell. Naturally, similar
keywords TwoProbe::RightElectrode::... exist for the right electrode.
To check that the alignment is made properly, the structure should be inspected visually.
When ATK is called with the command line option -c, it will stop after generating the
two-probe geometry. The xyz file can then be viewed using e.g. VIRTUAL NANOLAB.
$ atk li-ideal.atk -o li-ideal-check.nc -c --xyz --xyz-na 3 --xyz-nb 3
The --xyz-na and --xyz-nb options repeat the structure the specified numbers of
times in the x and y directions, and in this way it can be verified that the computational
cell is setup properly. The geometry is output in the file li-ideal.xyz.
When the atoms in the electrodes are aligned relatively to the atoms in the central
region, the two-probe cell vectors can be determined from the boundaries of the elec-
trode cells.
In the current version of ATK, the electrostatic problem is by default calculated by an

FFT routine, and the entire system must form a periodic structure. This means that the
first layer in left electrode must be similar to the first layer in the right electrode. It
is however possible to have different electrodes. In that case, one must include the
keyword TwoProbe::UseMultigridForElectrostatics in the input file, and as the
name indicates the Poisson equation is then solved by a multigrid method instead.
5.1 Geometry of an Ideal Lithium BCC Lattice 45
5.1.3 MISWRAPPED ELECTRODE ATOMS
A frequent setup mistake with the two-probe geometry is miswrapped atoms. Mis-
wrapping happens when the electrode atoms are aligned such that some atoms end
up outside the supercell. An example of miswrapping is illustrated in Figure 5.2. The
example is generated by adding an additional atom
Li 0.00 0.0 10.53
to the atom list in li-ideal.
5.1.4 TRANSMISSION SPECTRUM
To perform a two-probe calculation of li-ideal.atk, type the command
$ atk li-ideal.atk li-ideal.out
After some waiting, the self-consistent solution will be obtained and stored in the
NETCDF file li-ideal.nc. From this file various information about the system can
be obtained. To calculate the transmission spectrum generate the file trans.atk with
the content
and run the command
$ atk trans.atk -a li-ideal.nc trans.out
Inspection of the file trans.out shows ideal transmission of the lithium bands as ex-
pected.
Note, that small negative numbers can appear in the transmission spectrum, due to
round-off or other numerical errors in the integration in a complex plane (cf. also be-
low). The absolute values of these negative numbers can be reduced by reducing the
value of Analysis::TransmissionSpectrum::Eta, cf. page 90.
Figure 5.2: Miswrapped two-probe geometry of a Li (100) crystal. The right electrode is aligned
such that some atoms are outside the computational cell.
Figure 5.3: The contour used for the complex energy integral.
5.1.5 ENERGY INTEGRALS
The two-probe calculations are based on non-equilibrium Green’s function (NEGF)

techniques, where the electron density is resolved in energy space. To obtain the total
density from this spectral density, we need to perform an integral along the energy axis.
Part of the integral is performed on a complex contour, as illustrated in Figure 5.3. A
detailed description of the contour can also be found in Ref. 1.
In the complex plane, the functions are smooth and only few points are needed for
very accurate integration. The default parameters are conservative and will usually
be sufficient. However, if the value of TwoProbe::Contour::DeltaEMin is too low,
there may be bound states below the contour and this will result in charge missing
in the integration. This will typically manifest itself in a rapidly diverging value of "q"
during the SCF cycle (as reported in the log file), and what generally happens is that
the calculation converges (easily), but not to the correct result but rather to a situation
where there are no electrons in the central region any more.
If this is the case, it is necessary to increase TwoProbe::Contour::DeltaEMin. For an

estimate of the required value, one can perform a molecular calculation for the isolated
atoms to obtain the eigenenergies of the valence states. Subsequently, check that the
contour starts below these states by increasing the value of TwoProbe::Contour::
DeltaEMin accordingly. If the new value is much larger than the default, it is generally
a good idea to also increase TwoProbe::Contour::NumPoints in order not to lose
accuracy in the contour integration.
5.2 BENZENE RING COUPLED TO LITHIUM ELECTRODES
In this example, we will investigate the conductance of a benzene ring coupled with
two lithium fcc (111) electrodes. This system is similar to that considered in Ref. 12,
where gold electrodes were considered. That calculation takes much longer time, so
5.2 Benzene Ring Coupled to Lithium Electrodes 47
TwoProbe::Contour::DeltaEMin real EnergyUnit

The distance from the Fermi level to the energy where the complex contour starts.
Recommended value is twice the binding energy of the lowest eigenstate. If there are
electronic states below the complex contour, charge will be missing in the integration.
TwoProbe::Contour::RealAxis::PointDensity real EnergyUnit
The spacing between the points in the close-to-real axis part of the contour in the
voltage bias window.
TwoProbe::Contour::RealAxis::Eta real EnergyUnit
Small complex number added to the denominator in the real axis integration. A corre-
sponding parameter defined by the keyword (Analysis::TransmissionSpectrum::
Eta) is used for calculating the transmission spectrum.
Table 5.1: Important parameters for controlling the complex contour integration. For units and
default values, refer to the keyword reference section on pp. 86– 104.
for this example we choose lithium, which has only a single s band; the behavior of
the system still closely resembles gold, but is computationally cheaper. Figure 5.4
illustrates the system.
5.2.1 ELECTRODES
The lithium electrodes form a 3 3 structure, with 9 atoms in each layer. The left
electrode contains 3 layers with the stacking (BCA), while the right electrode has the
stacking (CBA). The stacking of the electrodes is chosen to match the central DTB (di-
thiol benzene) region which has the structure (BC–DTB–CBA). The complete two-probe
structure has the geometry (BCA)BC–DTB–CBA(CBA).
The complete two-probe structure is not perfectly periodic, but has a stacking fault at
the interface. This can be seen by repeating the structure (BCA)BC–DTB–CBA(CBA)|(BCA)BC–
DTB–CBA(CBA), which has a stacking fault (BA|BC) at the interface (the ideal interface
is BA|CB). Such small disturbances of the periodicity are allowed, and do not affect the
calculation.
In order to setup the electrode files, we will use the ATK keyword Electrode::NumRep.
Figure 5.4: A benzene molecule between two Li (111) surfaces. Note that some of the atoms in
the benzene ring have been periodically wrapped. Wrapping of atoms in the central region does
not affect the calculation, since periodic boundary conditions apply in the transverse directions.
With this keyword, the file for the left electrode (li111-bca.atk) becomes:
Method::Type NumOrb
Electrode::NumRep::A 3
Electrode::NumRep::B 3
UnitCell::LatticeConstant 5.4537 Bohr
0.500000 -0.866025 0.0000000
0.500000 0.866025 0.0000000
0.000000 0.000000 2.4494897
AtomList::Format Fractional
Li 0.3333333 0.66666666 0.0000000
Li 0.6666666 0.33333333 0.3333333
Li 0.0000000 0.00000000 0.6666666
The Electrode::NumRep parameter constructs a super cell consisting of the specified

unit cell repeated Electrode::NumRep times along the specified axes. Hence, in the
above example, we construct a supercell with 3 3 = 9 atoms per layer. The benefit
of setting up the system, which will contain 27 atoms in the end, through the keywords
Electrode::NumRep instead of typing in the coordinates of the 27 atoms, is that the
supercell calculation can be performed with 3 atoms and 9 k-points in the A/B-plane,
instead of with 27 atoms and 1 k-point. The results would be the same, as one can
convince oneself of by applying Bloch’s theorem and by folding the supercell Brillouin
zone accordingly, but the 27-atom case would take much longer since the calculation
time scales with the number of atoms to the third power, but only linearly with the
number of k-points. Moreover, when running in a parallel computing environment,
the program is designed to parallelize the k-point sampling, and so in this case the
27-atom calculation would, with access to 9 parallel nodes, take the same time as a
3-atom calculation.
ATK will automatically adjust the k-point grid to account for the zone folding, so
the k-point sampling that is specified in the electrode input file should be the de-
sired sampling for the complete system, i.e. ultimately the composite two-probe sys-
tem. To be specific, in the present example the final two-probe calculation will be
performed with a single k-point in the A/B plane, since we use the default values
Electrode::NumKPoints::A/B (note that the sampling in the C-direction is a com-
pletely different matter), while the actual electrode calculation will be performed with
9 points.
The use of the keywords Electrode::NumRep is limited to the electrode input files. In
5.2 Benzene Ring Coupled to Lithium Electrodes 49
the central region, it is necessary to explicitly list the coordinates of all atoms. Thus, if
we for example wish to include 3 surface layers of Li in the central region, we need to
give the coordinates of all 27 surface atoms in the two-probe input file.
In the right electrode file li111-cba.atk, the electrode layers are shifted to reflect the
stacking sequence CBA, i.e. it contains the following atom list
Li 0.3333333 0.66666666 0.0000000
Li 0.6666666 0.33333333 0.3333333
Li 0.0000000 0.00000000 0.6666666
All other lines are identical to li111-bca.atk.
5.2.2 CENTRAL REGION
To glue together the central region with the electrodes, we need to specify the equiva-
lent electrode and central atoms. By default, ATK will use the first and the last atom in
the central region atom list to match the first and last atom of the left and the right elec-
trodes, respectively. In the following example, we ensure this by placing the matching
atoms as the first and last atom in the central atom list:
Method::Type NumOrb
TwoProbe::LeftElectrode::ATKFile li111-bca.atk
TwoProbe::RightElectrode::ATKFile li111-cba.atk
#(BCA)BC-DTB-CBA(CBA)
Li 1.44298600 0.83310800 5.89096100
..
.
Li 0.00000000 0.00000000 23.00722200

The central atom list is too long to include here; it can be found in the file li-dtb-li+
0V.atk in the examples/twoprobe directory.
1.2
0V
1V
1
0.8
0.6
0.4
0.2
0
−4 −2 0 2 4
Energy (eV)
Figure 5.5: The transmission coefficient of the Li–DTB–Li system at 0 V and 1 V bias.
As opposed to the previous example, we do not make independent runs of the elec-
trodes. Instead, we use the keywords TwoProbe::LeftElectrode::ATKfile and TwoProbe::
RightElectrode::ATKfile. In this way, the electrode calculations are made implicitly
when performing the two-probe calculation.
Calculate the transmission spectrum with the command
$ atk li-dtb-li+0V.atk li-dtb-li+0V.out
The calculation will take some time. The transmission spectrum, extracted from the
output file, is illustrated in Figure 5.5.
5.2.3 FINITE BIAS
To perform a calculation at finite bias, we simply shift the chemical potential of the left
electrode relative to the right electrode. To do this, add the command
TwoProbe::LeftElectrode::Voltage 1 eV
and run the new file:
$ atk li-dtb-li+1V.atk li-dtb-li+1V.out -i li-dtb-li+0V.nc
With the -i option we use the zero bias density as a initial guess for the self-consistent
calculation to speed up the calculation.
5.3 Heterogeneous Electrodes 51
The result is shown in Figure 5.5. Due to the applied bias, the band edge of lithium has
moved closer to the HOMO level of the DTB, and part of the transmission is cut off at
low energies.
5.3 HETEROGENEOUS ELECTRODES
A new and unique feature in ATK is the possibility to have different electrode materials.
This can be used to study systems such as a metal–nanotube contact, or an interface
between two different metals.
There are a couple of separate aspects to keep in mind when setting up a heterogeneous
problem, compared to all other two-probe system so far considered in this tutorial,
where the electrodes were always the same (except for some differences in the stacking
sequences).
First of all, we mentioned earlier in this chapter that there is a requirement in ATK
that the electrodes must be periodic in all directions, including the transport or C
direction. The underlying reason for this requirement is that the Poisson equation,
which is an integral part of the self-consistent procedure, is most efficiently solved by a
Fourier-transform method. There is, however, another method for solving this equation
which does allow for the electrodes to be different. This method is called the multigrid
method, and it has now been implemented in ATK. By setting the keyword
TwoProbe::UseMultigridForElectrostatics
to True (by default it is False), the requirement of a periodic electrode cell in the C
direction is eliminated. The multigrid method is usually quite a bit slower than the
Fourier transform method, which is why it is not used by default for electrodes which
indeed are periodic.
The procedure of solving the Poisson equation is of course not restricted to the trans-
port direction. In the directions perpendicular to the interface plane, the periodicity
requirement still holds, which is very important to keep in mind. It means, for instance,
that if we want to consider an interface between two different metals, these must have
the same lattice constant. Since it is highly unlikely to find two such materials, one has
to resort to various tricks to circumvent this problem.
One way is it stretch one of the materials slightly (i.e. introduce strained materials),
an option which is viable if the difference in the lattice constants is small (such as for
instance between Au and Al). Another option is to try to come closer to commensurate
cells by making surface supercells, i.e. repeat the two surface unit cells a different
number of times, in the hope that, say, 3 cells of one material and 4 cells of the other
become almost the same size.
Figure 5.6: The geometry of the heterogeneous [111] bcc Li–C linear chain two-probe system,
visualized in VIRTUAL NANOLAB. Note that because we use the keywords NumRep to set up the
Li electrode, only the primitive electrode unit cell is shown. In the central region, however, we
have full layers of Li atoms.
Note that it is not sufficient to just come close; the two electrode cells must have exactly
same size in the A and B directions. Thus, usually one additionally has to include some
strain to get a perfect match.
One type of system where this latter requirement is trivial to fulfill is if one electrode
(or both) is a quasi one-dimensional system, such as a chain of atoms or a nanotube. In
this case the unit cell in the A and B directions can be chosen arbitrarily (as long as it is
large enough to avoid interactions between the repeated copies of the central region).
To demonstrate this, we shall here consider an example of a contact between a lithium

surface and a atomic carbon wire. This is admittedly a somewhat artificial system;
a more realistic example would instead be, say, a carbon nanotube and a gold or
palladium surface. The principles involved in setting up and calculating such a system
are however precisely the same as for our toy model. The only real difference is that
our model can easily be calculated on a standard PC within an hour, whereas a more
realistic system probably requires the use of a cluster (to take advantage of the parallel
power in ATK, in particular for k-point sampling), and uses much more memory. Still,
carbon wires are actually a topic of current research, at least when embedded inside
carbon nanotubes, so the system perhaps has some relevance anyway!
First, let us consider the left electrode, which is a Li bcc 33 [111] surface. The proper
way to set such a system up is by using the keywords Electrode::NumRep; that way we
only need to specify the three atoms in the primitive surface unit cell (li111-3x3.atk):
Electrode::NumKPoints::A 4
Electrode::NumKPoints::B 4
Electrode::NumRep::A 3
Electrode::NumRep::B 3
UnitCell::LatticeConstant 1. Angstrom
4.9629 0 0
-2.48145 4.298 0
0 0 6.07829
Li 2.48145 1.43267 0.506524
Li 0 0 1.51957
Li 0 2.86533 2.53262
Li 2.48145 1.43267 3.54567
Li 0 0 4.55871
Li 0 2.86533 5.57176
Three layers of Li are required to have a periodic cell for the [111] cleaving, and these
atoms are also sufficient to have a long enough electrode cell. In addition to run-
ning much faster, the NumRep calculations uses far less memory, which can be a real
concern on workstations. We use a small basis set, again to make the calculations
manageable and fast, but try to make the k-point sampling a bit more realistic by not
just using the default 11 sampling in the interface plane.
The right electrode is an atomic carbon wire. Precisely as we did for the Li chain in
Section 2.2 we must determine the spacing between the carbon atoms by minimizing
the total energy. We can basically use the same script as we did then (and for that
reason, this is not included in the scrips shipped with ATK); the result is that the minimal
energy is obtained with a distance of 1.34 Å between the C atoms (using a SZ basis
set).
In order to have a long enough electrode cell it is sufficient to include two repetitions,
at least to ensure that all first-order interactions are accounted for (cchain.atk):

14.8887 0. 0.
-7.44435 12.894 0.
0. 0. 2.68
C 2.48145 7.16337 0
C 2.48145 7.16337 1.34
Pay attention to the unit cells in the two electrodes file, and note that they are identical
in the A and B directions, once the 33 repetition is applied, as required.
Note also that we are forced to use the same k-point sampling in both electrodes,
despite the fact that we could have settled for only a single k-point in the A and B
directions for the carbon chain. This is because the sampling in the central region
is defined to be the same as in the electrodes (which is why there are no keywords
TwoProbe::NumKPoints), and hence the electrodes must have the same sampling.
Running the electrodes,
$ atk li111-3x3.atk li111-3x3.out

$ atk cchain.atk cchain.out
should be rather fast.
In the central region we put three layers of Li and 8 carbon atoms. This is most likely
enough for proper screening on both sides (we will soon verify this by inspecting the
results); linear chains are generally quite poor at screening surfaces (a pure geometrical
consequence) so one should expect to use a couple of extra layers on the chain side,
which fortunately not is particularly costly as each layer contains only a single atom.
We consistently use the smallest basis set (SZ) and a reduced mesh cut-off to speed up
the calculations.
Thus in total there are 35 atoms in the central region (li-cchain.atk):
NumOrb::MeshCutOff 100 Ry
TwoProbe::LeftElectrode::NetCDFFile li111-3x3.nc
TwoProbe::RightElectrode::NetCDFFile cchain.nc
TwoProbe::Contour::DeltaEMin 4. Ry
TwoProbe::Contour::Circle::NumPoints 50
SCF::TwoProbe::InitialBulkRun T
SCF::DiagonalMixingParameter 0.2
SCF::HistorySteps 20
TwoProbe::UseMultigridForElectrostatics T
Li 0.0000 0.0000 6.5848
Li -2.4815 -1.4327 7.5979
..
.
C 0 5.73067 16.26
C 0 5.73067 17.60
C 0 5.73067 18.94
C 0 5.73067 20.28
As for the other parameters used in the two-probe input file, an initial bulk run is of-
ten also useful for heterogeneous systems. We use a lot of history steps, otherwise it
will be difficult to reach convergence, and for safety we have increase the values of
TwoProbe::Contour::DeltaEMin and to avoid the risk of losing charges in the com-
plex contour integration.
We shall run the two-probe system using the --verbose option to get a bit more infor-
mation about the convergence progress:
$ atk li-cchain.atk li-cchain.out --verbose
The run will take about an hour. Considering that the total amount of k-points is 400, it
is clear that this system (both electrode and two-probe calculations) is extremely well
suited to calculate in parallel, and the performance should scale almost linearly with
the number of processors.
Inspecting the output file, we observe that the Mulliken populations change relatively
little during the second self-consistent cycle, indicating that the convergence is stable,
and that the initial bulk run did indeed provide a good initial guess for the calculation
with open boundary conditions. This is the reason why we can afford such a large
mixing parameter (we set SCF::DiagonalMixingParameter to 0.2 in the two-probe
run). We also see that the population of the C atoms quickly drops back to 4 as we
move away from the surface, but still
Figure 5.7: The left part of the figure shows a contour plot (cutting through the C atoms) of the
so-called electrostatic difference potential. This potential, which is defined as the solution to
the Poisson equation for an electron density corresponding to the difference between the atomic
densities and the self-consistent density, is relatively constant far away from the contact region,
both deep in the Li surface and in the C chain. The right plot show the effective single-electron
potential. Obviously the carbon chain extends to the left. The scale is different in the two plots,
but the colors work the same way.
Pay attention to the forces in the final result. The C atom closest to the surface expe-
riences an almost vanishing force in the z-direction (this was, in fact, the criterion for
choosing its distance to the surface), but the other atoms in the C chain are subject to
(large) forces that alternate in direction. Thus it is clear that in the presence of the Li
surface one must reconsider the position of all C atoms in the central region, and a
proper relaxation should be carried out. (There are also forces on the Li atoms in the
central region, but one can expect that the influence on the surface from the chain is
much smaller than the opposite effect.)
Let us look further at the Mulliken populations (some irrelevant lines are cut out).
%block Results::MullikenPopulations
0 Li q = 1.03155 [ s: 1.032 ]
1 Li q = 1.03915 [ s: 1.039 ]
2 Li q = 0.85846 [ s: 0.858 ]
3 Li q = 1.03078 [ s: 1.031 ]
4 Li q = 1.03918 [ s: 1.039 ]
..
.
13 Li q = 1.03556 [ s: 1.036 ]
14 Li q = 0.63550 [ s: 0.636 ]
15 Li q = 1.03078 [ s: 1.031 ]
16 Li q = 1.03616 [ s: 1.036 ]
17 Li q = 0.85840 [ s: 0.858 ]
..
.
27 C q = 4.27496 [ s: 1.510, y: 0.811, z: 1.144, x: 0.810 ]

28 C q = 4.18817 [ s: 0.961, y: 1.079, z: 1.068, x: 1.080 ]
29 C q = 4.07750 [ s: 0.942, y: 1.023, z: 1.089, x: 1.024 ]
30 C q = 4.07890 [ s: 0.919, y: 1.042, z: 1.076, x: 1.042 ]
31 C q = 4.05772 [ s: 0.922, y: 1.027, z: 1.081, x: 1.028 ]
32 C q = 4.04744 [ s: 0.920, y: 1.023, z: 1.081, x: 1.024 ]
33 C q = 4.03796 [ s: 0.918, y: 1.019, z: 1.081, x: 1.020 ]
34 C q = 4.05127 [ s: 0.927, y: 1.016, z: 1.092, x: 1.016 ]
# ----------------------------------------------------------------
# Total charge = 59.0256
# ----------------------------------------------------------------
%endblock Results::MullikenPopulations
Note how the bonding between the surface and chain is manifested in charge transfer
from the Li atom number 14 (the contact point) into the C chain. An isolated Li atom
has a population of 1 and a C atom 4, but in the converged two-probe system the
contact Li atom only has a population of 0.63, while the two closest C atoms in the
chain harbor about 4.2–4.3 electrons. One can also observe how the C atom closest
to the surface has its p-shell polarized along the z-direction, while all other atoms are
more symmetrically populated. Finally, the total Mulliken population approaches the
bulk value 4 as we move away from the surface, and the most distance Li atoms also
have populations close to unity, showing that the screening is sufficient.
The charge transfer between the Li and the C chain is not particularly large, and the
system has a very low conductance. One might imagine that the conductance would
be larger if a different contact point was used (such as an interstitial surface point,
where the chain could penetrate deeper). By moving the chain around one could
map out the conductance as a function of the contact point. The conductance would
most likely also change considerably if we were to relax the positions of the C atoms
more carefully, as mentioned above. All this is however well beyond the scope of this
tutorial, and so we leave it up to the interested reader to carry out this analysis!
C H A P T E R
6
SPIN-POLARIZATION AND k-POINT
SAMPLING OF TWO-PROBE
SYSTEMS
In this chapter, we focus on some of the more advanced features in the two-probe
calculations. The examples will illustrate how to treat spin polarization for two-probe
systems and the use of k-point sampling. Other issues that will be addressed is how to
achieve convergence for difficult cases and examples of how to extract data from the
NETCDF file.
The input files for these calculations are located in examples/twoprobe-spin in the
installation directory.
6.1 TRANSMISSION THROUGH A FE–MGO–FE INTERFACE
The example in this chapter illustrates the transmission through a Fe–MgO–Fe interface.
This is a system which currently recieves a lot of attention in the field of spintronics [3,
14], because it turns out to be an efficient spin filter, i.e. it only allows the majority spin
electrons of a current to pass through the interface.
In this chapter, we investigate a realistic application of ATK, and the calculations are
therefore rather time consuming. It takes about 24 hours to complete the calculations
in this chapter on a modern PC. Users with the possibility to run ATK in parallel should
however experience a substantial speedup on these examples.
59
60 Spin-Polarization and k-point Sampling of Two-Probe Systems
Figure 6.1: Geometry of the electrode and central region of the Fe-MgO-Fe system. Fe atoms are
turquise, O atoms red, and Mg atoms violet.
The system we will investigate is illustrated in Figure 6.1. The system has previously
been studied numericaly by Wortmann et. al. [15].
6.1.1 FE ELECTRODES
The Fe (100) electrodes are defined by the following input file (fe-100.atk)
Electrode::ElectronTemperature 0.1 eV
Method::Type NumOrb
NumOrb::BasisSet::Size SZP
NumOrb::XCFunctional LSDA-PZ
SCF::Tolerance 1.e-5
UnitCell::LatticeConstant 1. Bohr
5.42 0.00 0.00
0.00 5.42 0.00
0.00 0.00 16.26
AtomList::Format Bohr
Fe 1.355 1.355 1.355
6.1 Transmission Through a Fe–MgO–Fe Interface 61
Fe 4.065 4.065 4.065

Fe 1.355 1.355 6.775
Fe 4.065 4.065 9.485
Fe 1.355 1.355 12.195
Fe 4.065 4.065 14.905
By setting the exchange-correlation potential to LSDA-PZ the system is modelled as a

spin polarized system. The initial spin is by default chosen as the maximum polariza-
tion of each atom, and with the majority spin pointing up. The total initial spin can also
be specified with the keyword SCF::InitialSpinPolarization. For instance, setting
SCF::InitialSpinPolarization -0.1 fractional
reverses the sign of the majority spin.
Note, that we have increased the electron temperature in the electrode. This is a
computational trick that will help us in getting a converged result.
To calculate the effective potential of the electrode region we use the command.
$ atk fe-100.atk fe-100.out
Because the electrode region of the left and the right electrode is the same, we can use
the same electrode file for both.
6.1.2 FE-MGO-FE GEOMETRY
We next set up the central region with the commands (fe-mgo-fe.atk)
Method::Type NumOrb
NumOrb::BasisSet::Size SZP
NumOrb::XCFunctional LSDA-PZ
AtomList::Format Bohr
Fe 1.355 1.355 1.355
Fe 4.065 4.065 4.065
Fe 1.355 1.355 6.775
Fe 4.065 4.065 9.485
Fe 1.355 1.355 12.035
Mg 4.065 4.065 16.205
O 1.355 1.355 16.205

Mg 1.355 1.355 20.355
O 4.065 4.065 20.355
Mg 4.065 4.065 24.505
O 1.355 1.355 24.505
Fe 1.355 1.355 28.675
Fe 4.065 4.065 31.225
Fe 1.355 1.355 33.935
Fe 4.065 4.065 36.645
TwoProbe::LeftElectrode::NetCDFFile ./fe-100.nc
TwoProbe::RightElectrode::NetCDFFile ./fe-100.nc
Note the use of the keyword SCF::HistorySteps. With this keyword, we improves the
performance of the SCF routine, at the expense of slightly more memory usage since
the algorithm now uses information from 20 previous steps to make a new step.
To perform the SCF calculation use the command
$ atk fe-mgo-fe.atk fe-mgo-fe.out
6.1.3 TRANSMISSION SPECTRUM
Upon inspection of the output file, fe-mgo-fe.out, we notice that there are two sec-
tions with the transmission spectrum, corresponding to the spin-up and spin-down
component of the electronic system. The transmission coefficients of the spin up and
down electrons are shown in Figure 6.2. Notice how the transmission of spin down
electrons are strongly suppressed at the Fermi-level. Hence, the calculation demon-
strates that the MgO interface indeed functions as a spin filter.
6.2 k-POINT SAMPLING IN TWO-PROBE SYSTEMS
The Fe–MgO–Fe interface system has periodic boundary conditions in the two direc-
tions parallel to the interface (the A and B directions). In these directions, the wave-
functions are Bloch waves and we must average over all the k-vectors when calculating
physical quantities. In our two-probe calculation, two different quantities depend on
the number of k-points used in the A and B directions. The first quantity is the the
density matrix, the other is the transmission coefficient. Usually, it will require more
k-points to accurately calculate the transmission coefficient, compared to the density
matrix, and in ATK, it is therefore possible to specify different numbers of k-points for
the two quantities.
6.2 k-point Sampling in Two-Probe Systems 63
100
10−1
10−2
10−3
10−4
10−5
Spin Up
−6
Spin Down
10
−2 −1 0 1 2
Energy (eV)
Figure 6.2: Transmission coefficients of the spin up and down electrons for the calculation de-
fined in the file fe-mgo-fe.atk.
The number of k-points used for the density matrix is defined in the electrode input
file. To perform a two-probe calculation with a 6 6 k-point mesh in the A and B
directions, we must specify this mesh for the electrode file. We can do this by adding
to the file (fe-100.atk) following lines
and save it in fe-100-6x6.atk and subsequently execute the command
$ atk fe-100-6x6.atk fe-100-6x6.out
The result of the electrode calculation can now be compared with the calculation with-
out k-points. We see that all quantities, like the total energy, Mulliken populations (only
shown when running with the {-}{-}verbose option), etc have changed significantly.
This is a clear sign that an insufficient k-point sampling was used in the first calcula-
tion. To have a well converged k-point sampling, the change in the total energy/atom
should as minimum be less than 10 meV, and preferably in the order of 1 meV/atom.
If we test different k-point samplings for the electrode, we find that indeed the 6 6
produces well converged results. By increasing the mesh to a 8 8 mesh the total
energy only changes by 1 meV/atom. A calculation with a 4 4 mesh, on the other
hand, produces results that differ by 10 meV/atom.
Thus, the electrode calculations show that the 6 6 mesh is a good choice, and we
next include the 6 6 electrode files in the two-probe input file simply by copying the
file fe-mgo-fe.atk to fe-mgo-fe-6x6.atk and replace the lines
100
10−1
10−2
10−3
10−4
10−5
Spin Up
−6
Spin Down
10
−2 −1 0 1 2
Energy (eV)
Figure 6.3: Transmission coefficients at the point for the two-probe system in Figure 6.1. The
different curves are for calculations using 1 1, 2 2, 4 4, 6 6, and 8 8 k-points for the
density matrix calculation. The bold lines show the result of the 6 6 calculation, as defined in
the input file fe-mgo-fe-6x6.atk
TwoProbe::LeftElectrode::NetCDFFile ./fe-100.nc
TwoProbe::RightElectrode::NetCDFFile ./fe-100.nc
with
TwoProbe::LeftElectrode::NetCDFFile ./fe-100-6x6.nc
TwoProbe::RightElectrode::NetCDFFile ./fe-100-6x6.nc
The two-probe calculation will use the k-point mesh of the electrodes and we perform
the self-consistent two-probe calculation with the command
$ atk fe-mgo-fe-6x6.atk fe-mgo-fe-6x6.out
Note that we have not yet specified a k-point mesh for the transmission calculation,
and the transmission calculation therefore uses the default k-point mesh, which is a
1 1 mesh, i.e. the point.
In Figure 6.3 we show how the transmission spectrum at the point as function of the
number of k-points used for converging the two-probe calculation. The transmissions
for 6 6 and 8 8 are indistinguishable, indicating that the 6 6 mesh is sufficient.
6.3 k-point Sampling of the Transmission Spectrum 65
6.3 k-POINT SAMPLING OF THE TRANSMISSION SPECTRUM
We now have a converged density matrix and thereby a converged Hamiltonian of

the system. The next step is to obtain a k-point converged transmission spectrum. The
calculation of the transmission spectrum is a post processing step, and we will run atk in
the analysis mode using the converged Hamiltonian from the file fe-mgo-fe-6x6.nc.
The k-point sampling of the transmission coefficient is specified with the keywords
Analysis::TransmissionSpectrum::NumKPoints::A
Analysis::TransmissionSpectrum::NumKPoints::B
To calculate the transmission coefficient at the Fermi-level with a 20 20 k-point

mesh, we use the following input file (transmission.atk)
Analysis::TransmissionSpectrum::Calculate T
Analysis::TransmissionSpectrum::E0 0 eV
Analysis::TransmissionSpectrum::E1 +1 eV
Analysis::TransmissionSpectrum::Eta 1.e-5 eV
Analysis::TransmissionSpectrum::NumKPoints::A 20
Analysis::TransmissionSpectrum::NumKPoints::B 20
To perform the calculation, type the command
$ atk transmission.atk transmission.out -a fe-mgo-fe-6x6.nc
In the transmission.out file, we can read the k-point averaged transmission coeffi-
cient at the Fermi-level E = 0. We see that it is significantly different from the transmis-
sion coefficient at the point (compare the numbers with Figure 6.2).
The density matrix was well converged for a 6 6 k-point sampling, however, even
20 20 k-points are not enough to converge the transmission coefficient. This is illus-
trated in Figure 6.4, which shows the spin-up and spin-down transmission coefficient
calculated with different numbers of k-point. To obtain a converged result, more than
20 000 k-points are needed.
6.3.1 PLOTTING THE k-DEPENDENT TRANSMISSION
To understand why so many k-points are needed to converge the transmission coef-
ficient, it is instructive to plot the k-dependence of the transmission coefficient. All
the data used for calculating the k-point averaged transmission coefficient are stored
in NETCDF data file, transmission.nc. In future versions of VIRTUAL NANOLAB, it
10−1
10−2
10−3
Spin Up
−4
Spin Down
10
10 100 1000 10000 100000
Number of k−points
Figure 6.4: Transmission coefficients of the Fe-MgO-Fe system as function of the number of
k-points used for the transmission calculation.
will be possible to manipulate such data from the graphical user interface, however,
for now you will have to extract the data manually and plot them with your favorite
plotting program
Many programs like MatLab or Python has NETCDF support, and can be used rather
conveniently to extract data from the NETCDF file. However, in the following we will
use the command ncdump to extract the data. This is not an elegant solution but the
command is available on most platforms. For more information, see Appendix C.1.
To obtain a list of all the variables contained in the NETCDF file, type
$ ncdump -h transmission.nc
The variables we want to extract from the NETCDF file are TransmissionSpectrumTransmission\
_Up and TransmissionSpectrumTransmission\_Down. To extract the data use the
command
$ ncdump -v TransmissionSpectrumTransmission_Up,TransmissionSpectrumTransmission_Down transmission.nc
After the header part of the output, the spin up and down transmission spectra are
printed. There are 400 data elements in each variable. To get the values of the corre-
sponding k-points we extract the variable TransmissionSpectrumKpoints
$ ncdump -v TransmissionSpectrumKPoints transmission.nc
We get 400 values, corresponding to the kx and ky of 200 k-points. Note that a
(20 20) k-point mesh has 400 points, but ATK only calculates the transmission
6.3 k-point Sampling of the Transmission Spectrum 67
Figure 6.5: k-dependent transmission coefficient for majority (spin up) and minority spin (spin
down) at the Fermi energy for the Fe–MgO–Fe system. A (200 200) k-point grid covering the
entire Brillouin zone was used for the calculation.
at 200 of the points, as the transmission at the missing points can be obtained from
T( k) = T(k).
We may now pair the k-points with the transmission coefficients, and obtain the plot
shown in Figure 6.5. We see that there is a strong relation between the k-vector and the
transmission coefficient, and the behavior for the two different spin directions are rather
different. For the majority spin direction, the point gives the highest transmission,
whereas for the minority spin direction, only some special points in the Brillouin zone
can transmit the current.
Can such spin filtering be observed in real devices? To answer this question note that
we model a perfect interface, whereas in reality there will be some surface roughness.
The surface roughness destroys the conservation of the momentum parallel to the inter-
face (kx , ky ). To quantify the effect of surface roughness one must perform calculations
of systems that resemble the experimental situation. The effect of surface roughness
can easily be investigated with ATK by increasing the unit cell, and adding some fluc-
tuations in the interface. However, such calculations are beyond the scope of this
tutorial.
A P P E N D I X
A
ATK VS . TRANSIESTA-C
ATK constitutes a major upgrade over its predecessor TRANSIESTA-C in terms of per-
formance and features. Although there are some changes in the interface, these are
rather limited, and the users will hopefully agree that the changes serve to simplify the
input files and make the software easier to use. Thus, we are confident that users of
TRANSIESTA-C will feel right at home with this new generation of Atomistix’ electronic
structure software.
Here is a short list of new or updated features:
A general performance boost has been achieved due to a complete redesign of

the program, plus improvements in the basic algorithms. The program also uses
less memory now.
It is now possible to use different electrode materials in two-probe systems, by

setting the keyword TwoProbe::UseMultiGridForElectrostatics to T (the de-
fault is F, in which case the faster FFT method is used, as in TRANSIESTA-C.).
Spin polarization has been implemented.
Support for GGA has been added.
The mixing algorithms have been improved by adding the Pulay mixing method,
and the calculations now often converge faster and easier than before.
The code is parallelized (by using MPI, more specifically MPICH), which gives a
great performance boost in particular for two-probe systems.
There are no longer any archive (ark) files; instead, all the analysis results as well
as the data required to restart a calculation is saved in NETCDF files with a public
format, giving a much more portable and transparent interface to VNL and other
programs that can read NETCDF files (see Appendix C.1 for more information).
69
70 ATK vs. TranSIESTA-C
The program no longer uses ion files, instead all basis sets are generated on the
fly, which allows for much more flexibility for the users to optimize their basis
set.
There are a couple of important bug fixes. Earlier, negative transmission values
sometimes appeared due to the fact that the surface Green’s function did not con-
verge, but this has been fixed. Constraining atoms in relaxations now works prop-
erly, and the calculation of the partial density of states (PDOS) has been updated.
Finally, k-point sampling for the self-consistent cycle in two-probe calculations
plus k-point sampling for the transmission spectrum now works properly.
ATK ships with a completely new set of improved pseudopotentials that cover the
entire periodic table up to element 103 (Lr).
ATK uses the FLEXlm license manager, which allows proper license control with
floating licenses, running in parallel, etc. This has important consequences, since
all users who have not yet upgraded their licenses must do so in order to run ATK
2.0. Please see the Installation Guide for more detailed information.
It should be noted that ATK employs double precision floating point variables as op-
posed to the single precision used in in TRANSIESTA-C. Although this may appear to be
a minor modification, which should only effect the last decimals of the results, the sit-
uation is not always so simple. The self-consistent loop is a highly complex non-linear
problem, and as such very sensitive to small changes in the parameters. Therefore, ATK
may produce slightly different results compared to TRANSIESTA-C, even if all other pa-
rameters are identically the same (including, notably, the pseudopotential; see below).
However, the discrepancy should lie roughly within the tolerance of the self-consistent
cycle, meaning that the same differences in the final numbers could have been ob-
tained by just using a slightly different tolerance or mixing parameter.
There are other sources of potential differences in the results for systems calculated in
TRANSIESTA-C and ATK. For many elements, we have updated and improved the pseu-
dopotentials; it is of course still possible to run the calculations with the old potentials
(or, for that matter, any other pseudopotential which conforms to the format understood
by ATK). Finally, the default values of some keywords have been changed, so unless
you specify these keywords explicitly in the input file, a different value compared to
TRANSIESTA-C will be used. Most notably, the default mixing algorithm is now Pulay
by default, which should give faster convergence in almost all cases.
An important point to note is that ATK is not fully backward compatible with old
TRANSIESTA-C input files. Certain keywords have been changed, others removed (and,
of course, some have been added). Most of the deprecated keywords are related to
functions that now are performed automatically in the program, but these keywords
were seldom used in input files anyway. If such keywords are encountered in the in-
put files, the program will inform about it, and the offending keywords should just be
removed from the input files.
More importantly, there are keywords for which the syntax has changed, or the keyword
A.1 Changes to Keywords 71
itself has a new name. These changes have been made in order to improve the logical
structure of the input files and to make it easier for new users to learn the atk language.
If such outdated keywords are found in the input file, ATK will inform about the new
syntax/keyword name, to make it easier to migrate the users’ existing input files to the
new format. The two most commonly used keywords in this category are NumOrb::
BasisType which now should be NumOrb::BasisSet::Size and SCF::Beta which
has been renamed to SCF::DiagonalMixingParameter. For a complete list of changes
in the keyword syntax, see the sections further ahead in this chapter.
In addition to the changes in the input keywords, users will also notice that the output
file has a slightly new format, but it should still be simple enough to read off the results
of the calculation.
A.1 CHANGES TO KEYWORDS
Below is a list of all the keywords which have been deprecated, or for which the syntax
has changed.
A.1.1 KEYWORDS NO LONGER AVAILABLE
The following keywords cannot be used anymore:
Analysis::3DGrids::SaveEffectivePotential Done by default.
Analysis::3DGrids::SaveElectronDensity Done by default.
Analysis::3DGrids::SaveTo Only NetCDF output format is supported.
Analysis::3DLine::... Not available.
NumOrb::Accelerate::... Accelerated runs are no longer available.
NumOrb::BasisSet::IonFile The program no longer uses external ion files.
NumOrb::IonFileDirectory The program no longer uses external ion files.
NumOrb::RadialSampling::MaxRange Now calculated automatically.
SCF::ToleranceEbs Not available, instead use SCF::Tolerance.
SCF::ToleranceH Not available, instead use SCF::Tolerance.
SCF::ToleranceRho Not available, instead use SCF::Tolerance.
SCF::UseTmpFile Not available.

Simulation::AccumulateHistory Not available.
Molecule::TotalCharge Not available.
TwoProbe::Contour::UseTmpFile Not available.
A.1.2 KEYWORDS WITH NEW NAMES
The following keywords have new names, but otherwise the same syntax.
Bulk::FermiTemperature :
is now Bulk::ElectronTemperature.
Electrode::FermiTemperature :
is now Electrode::ElectronTemperature.
Molecule::FermiTemperature :
is now Molecule::ElectronTemperature.
NumOrb::BasisType :
is now NumOrb::BasisSet::Size.
NumOrb::MoleculePaddingFactor :
is now Molecule::ElectrostaticPaddingFactor.
NumOrb::PseudoFileDirectory :
is now NumOrb::BasisSet::PseudoPotentialDirectory.
SCF::Beta :
is now SCF::DiagonalMixingParameter.
SCF::NumStepsReset :
is now SCF::HistorySteps.
A.1.3 KEYWORDS WITH NEW NAMES AND SYNTAX
The following keywords have new names, and a slightly modified syntax. Please refer
to Appendix D.2.
SCF::MixH :
is now SCF::Coordinate with a new syntax.
TwoProbe::LeftElectrode::ArchiveFile :
is now TwoProbe::LeftElectrode::NetCDFfile or
TwoProbe::LeftElectrode::ATKfile.
A.1 Changes to Keywords 73
TwoProbe::RightElectrode::ArchiveFile :
is now TwoProbe::RightElectrode::NetCDFfile or
TwoProbe::RightElectrode::ATKfile.
Analysis::BandLine::... :
has a new syntax.
A P P E N D I X
B
BASIS SETS
B.1 BASIS SET PARAMETERS
ATK uses SIESTA type numerical orbital basis sets [11]. In Table B.1 we list the param-
eters used to define the SZ, DZ, DZP, DZDP basis sets. It is possible to change these
parameters in the input file and in this way tailor the basis set to a specific application.
In the following subsections we describe the meaning of these parameters and the basis
orbitals generation procedure.
B.1.1 SINGLE ZETA ORBITALS
The basis functions are found by solving the radial Schrödinger equation for the atom
with a confinement potential. The confinement potential is defined by the parameters,
V0 , rinn and rc through the equation
8
>
<0 if r
< rinn
Vconf (r) = V0 exp[ 1=(r rinn )]=(rc r) if rinn < r < rc
>
:1 if rc < r
Figure B.1 shows the confinement potential and the corresponding basis functions used
for constructing the ATK standard basis set for hydrogen.
75
76 Basis Sets
First zeta parameters

NumOrb::Basisset::EnergyShift (real EnergyUnit)
The energy shift defines the confinement radius, rc of the numerical orbital. The
confinement radius is taken at the position of the first node of the wave function with
energy Enl +EnergyShift, where Enl is the energy of the atomic eigenstate.
(Default value: 0.01 Ry)
NumOrb::Basisset::DeltaRinn (real)
Parameter used to determine the inner radius, rinn , of the soft confinement potential.
rinn is determined by rinn = DeltaRinn*rc .
(Default value: 0.8)

NumOrb::Basisset::V0 (real EnergyUnit)
Parameter for defining the softness of the confined potential. If V0 = 0, a hard wall
potential is used.
(Default value: 40. Ry)
NumOrb::Basisset::Charge (real)
Setting a net charge of the ion, when generating basis orbitals.
(Default value: 0.)
Second zeta parameters
NumOrb::Basisset::SplitNorm (real)
Defines the matching radius of an analytical orbital which splits the first zeta orbital
into a double zeta basis. The matching radius is determined by requiring that the
total norm of the split orbital should be a certain fraction of the norm of the first zeta
orbital.
Table B.1: Parameters for controlling the basis set generation.
B.1.2 DOUBLE ZETA ORBITALS
The double zeta orbital is constructed through a procedure reminiscent of the splitting
of a Gaussian basis set. The split orbital is obtained by constructing an analytical basis
orbital that matches the first zeta orbital smoothly at the radius rsplit . The functional
form used for the split orbital is
(l
r (al bl r2 ) if r < rlsplit
2l (r) =
1l (r) if r rl
split
The radius rsplit is determined by specifying the norm of the split orbital through the
keyword NumOrb::Basisset::SplitNorm.
B.1.3 POLARIZATION ORBITALS
The polarization orbitals have higher angular momentum than the valence orbitals. In
ATK such orbitals are generated by perturbing the single zeta orbitals by an electric
B.2 Pseudopotentials 77
Figure B.1: The lower part of the plot shows the l = 0 effective potential for hydrogen (dashed)
with the soft confinement potential (solid). The upper part shows the lowest occupied state of the
confined potential (solid line), and the atomic s-wavefunction is indicated by the dashed curve.
The dotted curve shows the solution with energy shift 0.01 Ry. The position of the first node of
this solution defines the position of rc .
field, and extracting the l + 1 component of the perturbed orbital.
B.2 PSEUDOPOTENTIALS
ATK uses norm-conserving pseudopotentials, and is shipped with a database of pseu-

dopotentials for the entire periodic table (cf. Appendix G). This database is constantly
updated and improved to provide accurate, general purpose pseudopotentials. How-
ever, some users may have their own pseudopotential databases which they would like
to use with ATK.
NumOrb::BasisSet::PseudoPotentialDirectory (string)
The string specifies the directory to look for pseudopotentials.
(Default value: $ATK_DATA_DIR/pseudopotentials/atk2.0)
NumOrb::BasisSet::PseudoPotentialFile::XX (string)
The string gives the name of the pseudopotential file for element Xx.
(Default value: Empty)
Table B.2: Parameters for controlling the pseudopotentials.
This is possible, since ATK can read several different pseudopotential formats. It reads
the unified pseudopotential format (upf) defined by the PWSCF consortium. From the
website http://www.pwscf.org it is possible to download tools that can convert sev-
78 Basis Sets
eral different pseudopotential formats into upf. ATK can also read the psf format used
by the SIESTA software.
In Table B.3 we list the pseudopotential formats compatible with ATK, and give refer-
ences for finding more information and sample pseudopotentials.
Extension HTML Reference

.cpmd http://www.cpmd.org/
.fhi http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP/
.psf http://www.uam.es/departamentos/ciencias/fismateriac/siesta/
.uspp http://www.physics.rutgers.edu/~dhv/uspp/
.upf http://www.pwscf.org/pseudo.htm
Table B.3: List of pseudopotential formats that can be used by ATK either directly (upf, psf) or
through conversion to the upf format. The web references point to the authors’ websites where
additional information on the format and sample pseudopotentials can be found.
A P P E N D I X
C
COMMAND LINE OPTIONS
The general syntax of a ATK call is

$ atk inputfile [outputfile] [OPTIONS]
where OPTIONS can be any combination of the following
--version
Print version number and exit.
-h, --help
Print available commands and exit.
--param "command1; command2; ..."

Equivalent to including the lines command1, command2, ... in the input file.
--force
Overwrite existing output files if any. Required if the files already exist, otherwise
the program will not proceed.
-o, --netcdf-file NetCDFfile.nc

Write data into file NetCDFfile.nc
-c, --check
Invoke checking mode. Do not perform simulation. Exit after checks are per-
formed. Note that this option produces a NETCDF file with the same name as the
input file (unless top option -o is used, of course).
-a, --analyze NetCDFfile.nc

Load NetCDFfile.nc and perform analysis.
79
80 Command Line Options
-i, --init-from NetCDFfile.nc

Initialize density matrix data from NetCDFfile.nc.
-v, --verbose
Output additional information about progress of simulation.
-q, --quiet
Output less information about progress of simulation.
--xyz
Output positions to .xyz file for visualization.
--xyz-file filename
Output xyz data to filename.
--xyz-na integer
Number of repetitions of the system along the A unit cell vector when generating
the .xyz file. The default value, used if this option is omitted, is 1.
--xyz-nb integer
Number of repetitions of the system along the B unit cell vector when generating
the .xyz file. The default value is 1.
--xyz-nc integer
Number of repetitions of the system along the C unit cell vector when generating
the .xyz file. The default value is 1.
C.1 INPUT AND OUTPUT FILES
The standard input and output files involved in a ATK run (the $ sign symbolizes the
command prompt)
$ atk filename.atk filename.out
are
filename.atk The input file with parameters for the calculation.
filename.out A text file containing the progress information and results of the calcu-
lation, along with all input data for the calculation, including parameters which
are set to default values. If the output file name is omitted, the output will be
written on the terminal, and can thus be piped to a file or inspected for progress
information. For more details about the contents of the output file, see Section E.
C.2 Environment Variables 81
filename.nc A binary file in the NETCDF format containing all data of the calcula-
tion. It can be used for restarting the calculation, as well as for analyzing the
simulation. Different types of data, like the electrostatic potential or the electron
density, are contained in this file and can be viewed with e.g. VIRTUAL NANOLAB.
To redirect data to another filename, use the command line option -o.
Note that the default behavior can be changed by using the command line options, in
order for instance to change the name of the output files. The output files are by default
(i.e. if the output file names are given without paths) in the directory from where the
program is run.
The NETCDF format is a platform independent binary format and therefore the file can
be copied between computers using different operating systems. For more informa-
tion about NETCDF files in general, it is recommended to have a look at the website
http://www.unidata.ucar.edu/software/netcdf/, where one also can download various
tools for managing NETCDF files, e.g. the program "ncdump" for extracting data (val-
ues) from NETCDF files. Moreover, NETCDF files can be read by other programs, e.g.
MATLAB and PYTHON, plus of course, as mentioned several times in this manual, VIR-
TUAL NANO L AB .
C.2 ENVIRONMENT VARIABLES
ATK_DATA_DIR This environment variable is used to locate the file containing the de-
fault ATK parameters, atk-defaults.atk. It is also used in the file atk-defaults.
atk itself to setup the path to the basis and pseudopotential files.
PATH It is recommended to include the ATK binary directory in your path.
For more details about how to set these variables up for different systems, please refer
to the Installation Guide.
82 Command Line Options
A P P E N D I X
D
PARAMETERS AND KEYWORDS
In this chapter all the keywords and parameters that ATK understands will be described
in detail. For each keyword, we will specify the syntax, list possible units and parame-
ters, and explain what these parameters mean.
The last section will clarify some points of the output (log) file which is generated during
a calculation.
D.1 SYNTAX
The parameters that control the setup of the system geometry and the calculations
performed with ATK are specified in an input file. Optionally, there is a possibility
to add keywords via the command line option --param. In either case, the general
specification of the parameters involves a keyword which ATK recognizes and relates
to a particular property, and a value of this property.
Some, but not all, values have units, and for those that do the unit is mandatory; there
are no default units (except for atomic coordinates, for which the unit is specified in a
different way; see below).
Each line in the input file looks like
keyword value [unit]
which is also the format for the --param commands. There are also parameters which
correspond to a table of values. These are set by
%block keyword
line 1 line 2 ...
83
84 Parameters and Keywords
%endblock keyword
Lines starting with the symbol # will not be read by the input parser, and can thus be
used for leaving out certain lines temporarily, or to introduce comments and documen-
tation of the calculational setup.
We shall in this Appendix consider first the values and units, and then go through the
keywords in alphabetical order. Note that ATK is case insensitive and the use of capital
letters in the documentation is only used in order to increase the readability.
In what follows, keywords are typeset in monospace, whereas unit specifications are
given in italic. When there are different options for specifying the value, these options
are separated by k. For example, AtomList::Format may be specified either by the
keywords Scaled or Fractional, or a LengthUnit such as Angstrom. Therefore, the
corresponding entry in the list below is given as
AtomList::Format LengthUnit k Scaled k Fractional
D.1.1 VALUES
A value may be another keyword (such as in the case of Method::Type which may be,
among other things, NumOrb) or a value may be a quantity of one of the following kinds:
real Floating point variable of the form 117.5 or 1.175e+2.
integer Integer value of the form 117
logical Logical of the form T or F (meaning true or false).
string A text string.
D.1.2 UNITS
All values describing a physical parameter must have specified a unit. In the current
version of ATK the basic output units are Å, eV and their combinations, e.g. eV/Å for
forces. In the output NETCDF file, however, the units are generally Bohr and Rydberg;
in each case, both in the log file and NETCDF files, all quantities are always given with
units, to avoid misunderstandings.
Valid units in the input files are
LengthUnit m, meter, nanometer, Bohr, a0, Ang, Angstrom

Unit of length (several alternative spellings are accepted, like "nanometres"); in
D.1 Syntax 85
SI units the keywords correspond to the following values:

Bohr, a0 0.52917710–10 m
nm, nanometer 10–9 m
Ang, Angstrom 10–10 m
InverseLengthUnit 1/m, 1/nanometer, 1/Bohr, 1/a0, 1/Ang, FracRecip
Unit of inverse length (see LengthUnit above for a more complete list). If the key-
word FracRecip is used, coordinates are given in units of the reciprocal lattice
vectors. In this case, the unit cell must be given in the input file, also for analysis
runs! In addition, the keyword System::Bulk must be included.
EnergyUnit eV, meV, Ha, Hartree, J, Joule, K, Kelvin, Ry, Rydberg

Unit of energy (some alternative spelling are also accepted, like "millielectron-
volts"); in SI units the keywords correspond to the following values:
eV 1.6021910–19 Joule
meV 1.60210–22 Joule
Ha, Hartree 4.35986341210–18 Joule
Ry, Rydberg 2.17993170610–18 Joule
K, Kelvin 1.380710–23 Joule
ForceUnit Ry/Bohr, Ha/Bohr, eVPerAng, eV/Ang, nN, nanoNewton
Unit of force (alternative spellings like "RyPerBohr" or "Hartree/Bohr" are also
accepted); in SI units the keywords correspond to the following values:
Ry/Bohr 41.19475536 Newton
Ha/Bohr 82.38951072 Newton
eV/Ang 1.6021910–9 Newton
nN, nanoNewton 10–9 Newton
TimeUnit s, seconds, AUT, fs
Unit of time; in SI units the keywords correspond to the following values:
AUT 2.41886587810–17 seconds
fs 10–15 seconds
SpinUnit HBar, fractional
Unit of spin; we use atomic units (~ or Hbar) or measure the total electron spin
in fractions of the total number of electrons (fractional). In atomic units the spin
of one electron is 12 . The degeneracy of the state is (2S+1), where S is the total
spin. Thus,
0.0 Hbar Singlet
0.5 Hbar Doublet
1.0 Hbar Triplet
1.5 Hbar Quartet
2.0 Hbar Quintet
In some cases, in particular in analysis mode, spin is gives as keywords Up/Down
to distinguish between e.g. the spin up or down molecular eigenvalues.
ConductanceUnit Siemens, G0
Unit of conductance; in SI units the keywords correspond to the following values:
G0 7.7480691710–5 Siemens
D.2 KEYWORDS
Keywords are divided into different groups or classes according to a hierarchial system,
like in a menu system with subentries. Each class of keywords have the same base-
name, and the basename signifies which part of the calculation the keyword controls.
The different levels in the hierarchy are separated by double colons ::. For exam-
ple, the keyword Bulk::NumKPoints::A belongs to the major class Bulk, which has a
subclass NumKPoints.
The major keyword classes are:
Analysis:: Controls the analysis of the simulation.
AtomList:: Controls how the atomic coordinates are defined in the input file.
Bulk:: Parameters for systems with bulk symmetry.
Electrode:: Parameters of the electrode region of a two-probe system.
Method:: Select which method to use for the electronic structure calculation.
Molecule:: Parameters for systems with molecular symmetry.
NumOrb:: Parameters of the DFT numerical orbital basis set.
Relaxation:: Parameters for the relaxation of the atomic coordinates.
SCF:: Parameters for the self-consistency loop.
Simulation:: Determine what type of simulation to perform.
System:: Determines the overall symmetry of the system.
TwoProbe:: Parameters for a two-probe system.
UnitCell:: Defines how the unit cell of a periodic system is defined in the input file.
As pointed out earlier in this manual, parameters not explicitly listed in the input atk
file will be given values from the defaults file (global or user-defined; cf. page 81).
There are furthermore certain required keywords for each system type, without which
it is not possible to perform the calculation. These keywords are rather few, and
quite logical. Naturally, for any calculation it is necessary to define the System::Type
and to give an Atomlist. In addition, it is necessary to specify a UnitCell for bulk
and electrodes, and finally TwoProbe::LeftElectrode::NetCDFfile or TwoProbe::
LeftElectrode::atkfile (and correspondingly for the right electrode) are required
for two-probe systems.
Below all the keywords in each class will be described in detail in terms of syntax. The
keywords are given in alphabetical order, and for each one is also listed the default
value used if the keyword is not explicitly given in the input file.
D.2 Keywords 87
D.2.1 ANALYSIS::
This class of keywords controls the analysis of the simulation performed after the sim-
ulation has been performed.
Analysis::BandLine::X::Start
Analysis::BandLine::X::End (real, real, real) InverseLengthUnit
Print the band structure along the line from Analysis::BandLine::X::Start to
Analysis::BandLine::X::End with the number of points defined by Analysis::
BandLine::X::NumPoints. X can be any letter in the alphabet and is used to
identify the bandlines when several lines are included in the calculation.
Analysis::BandLine::X::NumPoints integer
Number of points along the band structure line associated with the label X. If not
specified, then the value Analysis::BandLine::NumPoints is used for all band
lines.
Analysis::BandLine::NumPoints integer
Number of points along the band structure lines.
(Default value: 100)
Analysis::BlochState::X::Index integer
Specification of the band index of the Bloch state to be calculated. X can be any
letter in the alphabet and is used to identify the Bloch state when several states
are included in the calculation. The numbering of the bands starts from 0.
Analysis::BlochState::X::KPoint (real, real, real) InverseLengthUnit

Specification of the wave vector, i.e. the k-point, of the Bloch state associated
with the index X (see Analysis::BlochState::X::Index).
(Default value: -point)
Analysis::BlochState::X::Label string
String used to label the Bloch state associated with the index X (see Analysis::
BlochState::X::Index).
Analysis::BlochState::X::Spin Up k Down
Specification of the spin of the Bloch state associated with the index X (see
Analysis::BlochState::X::Index). Only relevant for spin-polarized calcula-
tions.
(Default value: Up)
Analysis::DOS::X::Energy real EnergyUnit

Specifies the energy, at which a spatially resolved density of states, n(E; r; k; ),
should be returned. This can in turn be used to make 3D plots. Currently, this
works only for two-probe systems. The spatially resolved DOS is defined by:
X
n(E; r; k; ) = j ;k; (r)j (E ";k; );

where E is the energy, specified by the keyword Analysis::DOS::X::Energy;

k is a point in reciprocal space, specified by the keyword Analysis::DOS::X::
KPoint; and is the spin specified by Analysis::DOS::X::Spin. The latter
keyword should only be specified for spin-polarized calculations. ;k; (r) is
an eigenstate with eigenenergy ";k; . X can be any letter in the alphabet and
is used to identify the spatially resolved DOS with this particular energy (and k-
point and spin).
Analysis::DOS::X::KPoint real, real, real InverseLengthUnit

Specifies a point in reciprocal space, where the spatially resolved DOS associated
with index X should be determined (see Analysis::DOS::X::Energy).
(Default value: -point)
Analysis::DOS::X::Label string
The string used to label the spatially resolved DOS associated with the index X
(see Analysis::DOS::X::Energy).
Analysis::DOS::X::Spin Up k Down
Specifies the spin of the spatially resolved DOS associated with index X (see
Analysis::DOS::X::Energy. Only relevant for spin-polarized calculations.
(Default value: Up)
Analysis::DOS::X::Type Total k Left k Right k NonEquilibrium

This keyword controls which part of the real space density of states (DOS) that
is saved. Total means the total DOS, Left and Right mean the DOS from
the left and right electrodes, respectively, and NonEquilibrium means the DOS
weighted with the Fermi occupation, i.e. the occupied eigenstates.
Analysis::MolecularOrbital::X::Index integer
Specification of the state index of the molecular eigenstate to be printed in the
analysis output file. X can be any letter in the alphabet and is used to identify
the molecular state, when several states are included in the calculation. The
numbering of the states starts from 0.
Analysis::MolecularOrbital::X::Label string
String used to label the molecular state associated with the index X (see Analysis::
MolecularOrbital::X::Index).
D.2 Keywords 89
Analysis::MolecularOrbital::X::Spin Up k Down
Specifies whether the orbital associated with index X should be a spin up or down
orbital. Only relevant for spin-polarized calculations.
(Default value: Up)
Analysis::MPSH::Calculate logical
The eigenvalues of the Molecular Projected Self-consistent Hamiltonian (MPSH)
are calculated if this parameter is set to T. The MPSH is defined as the Hamilto-
nian of a two-probe system excluding the electrode and surface atoms. The num-
ber of surface atoms is determined by the two parameters TwoProbe::LeftSurface::
NumberOfAtoms and TwoProbe::RightSurface::NumberOfAtoms.
(Default value: T)
Analysis::MPSH::EnergyWindow real EnergyUnit

If this parameter is positive (i.e. > 0), only eigenvalues in an energy window of
the specified size around the average electrode Fermi level are reported.
(Default value: 0. eV)
Analysis::MPSH::X::Index integer
Specification of the state index of the MPSH eigenstate to be printed in the anal-
ysis output file. X can be any letter in the alphabet and is used to identify the
molecular state when several states are included in the calculation. The number-
ing of the states starts from 0.
Analysis::MPSH::X::Label string
String used to label the MPSH state associated with the index X (see Analysis::
MPSH::X::Index).
Analysis::MPSH::X::Spin Up k Down
Specifies whether the MPSH eigenstate associated with index X should be a spin
up or down state. Only relevant for spin-polarized calculations.
(Default value: Up)
Analysis::PrintMolecularEigenvalues logical
Controls if the molecular eigenvalues are printed in the log file or not.
(Default value: T)
Analysis::PrintMullikenPopulation logical
Perform Mulliken population analysis. The results are only displayed if the com-
mand line option --verbose is used.
(Default value: T)
Analysis::TransmissionSpectrum::Calculate logical
Calculate the transmission spectrum.
(Default value: T)
Analysis::TransmissionSpectrum::E0 real EnergyUnit

Value of first energy point of the transmission spectrum.
(Default value: -2. eV)
Analysis::TransmissionSpectrum::E1 real EnergyUnit

Value of last energy point of the transmission spectrum.
(Default value: +2. eV)
Analysis::TransmissionSpectrum::Eta real EnergyUnit

Small complex number added to the energy when calculating the transmission
spectrum.
(Default value: 1.e-5 eV)
Analysis::TransmissionSpectrum::NumEigenChannels integer
The transmission coefficient can be decomposed into transmission eigenchannels
by diagonalizing the transmission matrix ty t. The number of channels reported are
controlled by the keyword Analysis::TransmissionSpectrum::NumEigenChannels,
and the eigenvalues are reported for the eigenchannels with the highest transmis-
sion.
(Default value: 3)
Analysis::TransmissionSpectrum::NumKPoints::A integer
Number of k-points used for the transmission spectrum in the A direction.
(Default value: 1)
Analysis::TransmissionSpectrum::NumKPoints::B integer
Number of k-points used for the transmission spectrum in the B direction.
(Default value: 1)
Analysis::TransmissionSpectrum::NumPoints integer
Number of energy points used for the transmission spectrum.
D.2.2 ATOMLIST::
This class controls how the atomic coordinates are read into ATK. The actual coordi-
nates are entered in a list, where each line defines an atom and its position, using four
entries:
%block keyword
atomname x y z
...
atomname x y z
%endblock keyword
where atomname is the symbol of the relevant chemical element (such as Au for a
gold atom), and keyword should be replaced by the respective keyword for the rel-
D.2 Keywords 91
evant system type: Bulk::AtomList for periodic systems, Molecule::AtomList for

molecules, Electrode::AtomList for the electrodes in a two-probe system, and fi-
nally TwoProbe::CentralAtomList for the scattering region in a two-probe system.
How the three real numbers x, y, and z on each line are interpreted is controlled by
the following two keywords:
AtomList::Format LengthUnit k Scaled k Fractional

If a LengthUnit is specified, then the coordinates x, y, and z are assumed to be
given in the specified unit of length (such as Angstrom). If the keyword Scaled
is used, the coordinates are assumed to be given in units of AtomList::Scale.
Finally, if the keyword Fractional is used, then the three real numbers x, y, and
z are assumed to specify the atomic positions r through the relation r = xa1 +
ya2 + za3 , where a1 , a2 , a3 are the three lattice vectors, specified by the relevant
UnitCell keyword.
(Default value: Ang)
AtomList::Scale real LengthUnit

Used to scale the atomic coordinates if AtomList::Format is set to Scaled.
D.2.3 BULK::
This class specifies parameters which are specific to systems with bulk symmetry.
Bulk::AtomList %block
The coordinates of the atoms in the basis, given in a list (c.f. AtomList on
page 90).
Bulk::ElectronTemperature real EnergyUnit

The electronic temperature of the bulk system.
(Default value: 300. K)
Bulk::NumKPoints::A integer
Number of evenly spaced k-points in the A direction. ATK uses the Monkhorst–
Pack method for choosing the mesh points. [6]
(Default value: 1)
Bulk::NumKPoints::B integer
Number of evenly spaced k-points in the B direction.
(Default value: 1)
Bulk::NumKPoints::C integer
Number of evenly spaced k-points in the C direction.
(Default value: 1)
Bulk::UnitCell %block
The three lattice vectors or unit cell vectors a1 , a2 and a3 , expressed in units
of UnitCell::LatticeConstant in Cartesian coordinates. For example, a bcc
lattice with 2 Å lattice constant is most conveniently specified by the following
input:
UnitCell::LatticeConstant 2. Angstrom
-0.5 0.5 0.5
0.5 -0.5 0.5
0.5 0.5 -0.5
This way of doing it allows the lattice constant to be changed with a single mod-
ification of the file, instead of having to recalculate all components of all unit
vectors. It is of course still possible to give the same unit cell as
-1 1 1
1 -1 1
1 1 -1
since the default value of UnitCell::LatticeConstant is 1 Å.

D.2.4 ELECTRODE::
This class defines the electrode region of a two-probe system. The effective potential in
the electrode region is obtained from a corresponding bulk system, and the electrode
parameters set up the connection with this corresponding system.
Electrode::AtomList %block
The coordinates of the atoms in the electrode, given in a list (c.f. AtomList on
page 90).
Electrode::ElectronTemperature real EnergyUnit

The electronic temperature of the electrode.
Electrode::NumKPoints::A integer
Number of evenly spaced k-points in the A direction. ATK uses the Monkhorst-
Pack method for choosing the mesh points [6]. The k-point sampling specified
for the electrode in the A and B direction is both used for the electrode calcula-
tion and the two-probe calculation. Since the same k-point sampling is used in
D.2 Keywords 93
the two calculations, usually a converged calculation can be obtained with few
points. Typically, a more dense sampling is needed for calculating the transmis-
sion function, see Analysis::TransmissionSpectrum::NumKPoints::A.
(Default value: 1)
Electrode::NumKPoints::B integer
Number of evenly spaced k-points in the B direction.
(Default value: 1)
Electrode::NumKPoints::C integer
Number of evenly spaced k-points in the C direction. This parameter is used
when calculating the self-consistent effective potential of the corresponding bulk
system. The two-probe system is infinite in the C-direction and for a consistent
electrode, it is important that its electronic structure accurately describes the
infinite system. For this reason a rather large value of Electrode::NumKPoints::
C must generally be chosen.
Electrode::NumRep::A integer
Number of replicas of the corresponding bulk system in the A direction needed to
match the A direction of the two-probe system. Note, that when using Electrode::
NumRep::A >1, the k-point sampling specified with Electrode::NumKPoints::A
is the k-point sampling for the two-probe system, i.e. the repeated electrode cell.
The corresponding k-point sampling of the non-repeated electrode cell is auto-
matically calculated by ATK.
(Default value: 1)
Electrode::NumRep::B integer
Number of replicas of the corresponding bulk system in the B direction needed
to match the B direction of the two-probe system.
(Default value: 1)
Electrode::UnitCell %block
Similar to Bulk::UnitCell.
D.2.5 METHOD::
Parameters related to the method used to calculate the total energy and forces.
Method::Type NumOrb The energy method used for the calculation. NumOrb corre-
sponds to using density functional theory (DFT) with the electronic wave func-
tions expanded in a basis set of atom-based numerical orbitals.
(Default value: NumOrb)
D.2.6 MOLECULE::
This class specifies parameters which are specific to systems with molecular symmetry.
Molecule::AtomList %block
The atomic coordinates for the molecule, given in a list (c.f. AtomList on page 90).
Molecule::ElectrostaticPaddingFactor real
A molecule is represented within a periodic cell with side lengths a, b and c.
The numerical orbitals have finite radius, and a0 , b0 , and c0 are defined as the
lengths where the molecule is not interacting with any of its periodic images in
each respective direction. The side length of the computational cell is defined by
a = (1+Molecule::ElectrostaticPaddingFactor) a0 ,
and equivalently for b and c.
Molecule::ElectronTemperature real EnergyUnit

The electronic temperature of the molecular system.
Molecule::TotalSpin real SpinUnit

By default the spin of the system is determined self-consistently. With the Molecule::
TotalSpin parameter it is possible to fix the spin of a molecule to a specific
value. The spin is given in atomic units (~), in which the spin of a single elec-
tron is 0.5. Note that a spin-polarized method (see NumOrb::XCFunctional) is
required when using non-zero spin values. We recommend to set the SCF::
InitialSpinPolarization to the same value as this parameter.
D.2.7 NUMORB::
This class specifies the parameters of the DFT numerical orbital basis set used to solve
the electronic structure problem.
NumOrb::BasisSet::Size::Xx SZ k DZ k SZP k DZP k DZDP

This keyword describes the size of the numerical orbital basis for element Xx,
where Xx is Li for lithium, for example. The corresponding keyword without the
::Xx extension gives a global value for all elements.
SZ A single basis orbital for each valence orbital.

DZ Two basis orbitals for each valence orbital.
SZP SZ and one basis orbital for the first unoccupied shell.
D.2 Keywords 95
DZP DZ and one basis orbital for the first unoccupied shell.
DZDP DZ and two basis orbitals for the first unoccupied shell.
The basis set is generated according to the valence shells listed on pp. 121–124.
If semicore states are included in the pseudopotential, then separate zeta basis
orbitals are generated for the semicore states.
(Default value: DZP)
NumOrb::BasisSet::Charge::Xx (real)
Defines the net charge of element Xx, when generating basis orbitals. The key-
word without the ::Xx extension gives a global value for all elements.
(Default value: 0.)
NumOrb::Basisset::DeltaRinn::Xx (real)
Determines the inner radius, rinn , of the soft confinement potential for element
Xx. The keyword without the ::Xx extension gives a global value for all elements.
rinn is determined by rinn = DeltaRinn*rc . Valid numbers are in the range
[0.0,1.0].
NumOrb::Basisset::EnergyShift::Xx (real EnergyUnit)

Determines the confinement radius rc of the numerical orbitals of element Xx.
The keyword without the ::Xx extension gives a global value for all elements.
The confinement radius is taken at the position of the first node of the wave
function with energy Enl +EnergyShift, where Enl is the energy of the atomic
eigenstate.
(Default value: 0.01 Ry)
NumOrb::Basisset::SplitNorm::Xx (real)
Defines the relative norm of the second zeta, relative to the first zeta for element
Xx. The keyword without the ::Xx extension gives a global value for all ele-
ments. The splitnorm is used to find the matching radius of an analytical orbital
which splits the first zeta orbital into a double zeta basis. The matching radius
is determined by specifying that the norm of the split orbital relative to the norm
of the first zeta orbital should have the value NumOrb::Basisset::SplitNorm.
Valid numbers are in the range [0.0,1.0].
NumOrb::Basisset::V0::Xx (real EnergyUnit)

Parameter for defining the softness of the confined potential for element Xx. The
keyword without the ::Xx extension gives a global value for all elements. If
V 0 = 0, a hard wall potential is used.
NumOrb::BasisSet::PseudoPotentialDirectory string
The string describes the directory in which to look for pseudopotential files.
(Default value: $ATK_DATA_DIR/pseudopotentials/atk2.0)
NumOrb::BasisSet::PseudoPotentialFile::XX (string)
The string gives the name of the pseudopotential file for element Xx.
NumOrb::MeshCutoff real EnergyUnit

This parameter controls the fineness of the grid used for evaluating real space
integrals and for solving the Poisson equation. The real space resolution x is
p
related to the mesh cutoff Ecut through the relation x = / Ecut , if Ecut is in
Rydberg and x in Bohr.
NumOrb::XCFunctional LDA-PZ k LSDA-PZ k GGA-PBE k SGGA-PBE

This parameter defines the approximation used for the exchange–correlation func-
tional.
LDA-PZ The local density approximation (LDA) with the Perdew–Zunger parametriza-
tion [8] of the correlation energy of a non spin-polarized homogenous elec-
tron gas calculated by Ceperly–Alder [4].
LSDA-PZ The local spin density approximation (LSDA) with the Perdew–Zunger
parametrization [8] of the correlation energy of a spin-polarized homoge-
neous electron gas calculated by Ceperly–Alder [4]. .
GGA-PBE The non spin-polarized generalized gradient approximation (GGA) in
the Perdew–Burke–Enzerhofen approximation. [7].
SGGA-PBE The spin-polarized generalized gradient approximation (GGA) in the
Perdew–Burke–Ernzerhof approximation. [7].
If the functional describes spin polarization, the spin of the system will be deter-
mined self-consistently. For molecular systems it is possible to fix the spin using
Molecule::TotalSpin.
(Default value: LDA-PZ)
NUMORB:: SPECIAL PARAMETERS
This set of parameters defines the radial representation of the numerical orbitals. We
recommend to use the default values of these parameters.
NumOrb::RadialSampling::DR real LengthUnit

The separation between the points used for the radial representation of the nu-
merical orbitals.
(Default value: 0.005 Bohr)
NumOrb::Tables::Cutoff real EnergyUnit

Two center integrals are done in Fourier space, and this parameter controls the
maximum Fourier component used for the integral.
D.2 Keywords 97
NumOrb::Tables::NumPoints integer
Controls the number of points used for the two-dimensional table, which contains
the two center integrals.
D.2.8 RELAXATION::
These keywords control the relaxation of the ionic coordinates.
Relaxation::ForceTolerance real forceUnit

The relaxation is converged when the magnitude of the largest force component
is lower than this value.
(Default value: 5.e-2 eV/Ang)
Relaxation::MaxSteps int
The relaxation makes at most MaxSteps iterations.
Relaxation::MaxDisplacement real LengthUnit

Sets the maximum allowed length of each step in the relaxation routine.
(Default value: 0.1 Bohr)
Relaxation::Type (CG)
Specifies the algorithm used in the minimization. At present, only the conjugate
gradient method is available.
(Default value: CG)
D.2.9 SCF::
This class of keywords controls the self-consistency loop.
SCF::Coordinate DensityMatrix k Hamiltonian k EffectivePotential

This keyword determines which physical quantity that is mixed between iterations
in the self-consistent cycle.
(Default value: Hamiltonian)
SCF::Algorithm Pulay k Broyden

This keyword controls which mixing algorithm to use. In most cases, in particular
for two-probe systems, Pulay mixing is more stable and faster than Broyden.
(Default value: Pulay)
SCF::Tolerance real
When the relative change of the density and the Hamiltonian matrix is less than
this value between two iterations, the self-consistent cycle is considered con-
verged.
(Default value: 1.e-4)
SCF::Criterion Strict
The criterion for when self-consistency is reached. The keyword Strict indicates
that all three convergence criteria (density matrix, Hamiltonian and energy of all
occupied states) must reach the required tolerance.
(Default value: Strict)
SCF::DiagonalMixingParameter real
This parameter defines the fraction of the new value of the SCF::Coordinate
which is mixed with the old value.
SCF::HistorySteps integer
Number of previous steps used in the mixing algorithm to generate the new state.
Linear mixing corresponds to setting SCF::HistorySteps to unity, but this is
highly unstable and the calculations usually never converge. In case of conver-
gence problems this parameter can be increased to include more of the previous
steps to find the optimal value of the new state. This requires more memory, but
increasing fileSCF::HistorySteps up to 20 or even a bit more should be safe on
most systems.
(Default value: 6)
SCF::InitialSpinPolarization real SpinUnit

Sets the initial total spin of the system to a specified value. The parameter only
has an effect for a spin polarized calculation (see NumOrb::XcFunctional).
If this parameter is not specified, the spin occupation of each orbital is set so
that the initial polarization is equal to the maximum atomic polarization of each
element. In formulas, this means that
o";i = min(Oi ; 1); (D.1)
o#;i = max(Oi 1; 0): (D.2)
where o";i is the occupation of the "up" spin channel in orbital i, and o#;i is the
occupation of the "down" channel in orbital i, and Oi = o";i + o#;i is the atomic
occupation of the shell. Since Oi can only be 0, 1 or 2, the two occupations oi
can take the values 0 (unoccupied orbital) and 1 (occupied orbital).
To be more concrete, we can consider a specific example. In its ground state,
the d-shell of Fe has the configuration 3d6 ["#]["]["]["]["], where the [: : :] boxes
symbolize the occupancies of the 5 degenerate d-states. Thus, o";3d 2 = 1,
o#;3d 2 = 1, o";3d 1 = 1, o#;3d 1 = 0, and so on.
With the keyword SCF::InitialSpinPolarization it is possible to choose an-
other initial polarization, which is applied to all atoms alike. The keyword defines
D.2 Keywords 99
the parameter , which re-scales the polarized state. The value of is calculated
from the initial spin-polarization,
X
= (SCF :: InitialSpinPolarization)= (o";i o#;i )
i
The initial spin polarization is now chosen as
o~";i = o";i + (1 )Oi =2 (D.3)
o~#;i = o#;i + (1 )Oi =2 (D.4)
Note that if > 1 then we set = 1, and if < –1 then we set = –1.
Again, an example illustrates the point. The initial spin polarization of the ground
state of Fe as given above is 4; the spin polarization is defined as
X
(o";i o#;i :
i
If instead we would like to have an initial polarization of 2, we set the key-

word SCF::InitialSpinPolarization to 1 Hbar. Note carefully that the spin
of one electron is ~=2! This gives = 1=2, and the orbital occupations become
o";3d 2 = 1, o#;3d 2 = 1, and o";i = 3=4, o#;i = 1=4 for all other d-orbitals.
It may seem weird to have a fractional occupation of the orbitals, but this is
not a problem; the occupations only gives the probability of finding the electron
somewhere, and this can of course be different from an integer (cf. the situation
with a finite electron temperature, where the HOMO and LUMO occupations
may be fractional as well).
SCF::MaxSteps integer
Maximum number of steps in the self-consistency loop.
SCF::TwoProbe::InitialBulkRun logical
If this parameter is set to true, a two-probe system will initially be treated as a
periodic (bulk) system, subject to periodic boundary conditions, and a full self-
consistent run will be performed on this system. The purpose of this is to generate
a better initial guess of the density matrix for the full two-probe SCF loop. If the
initial bulk run is left out, the initial density matrix is evaluated from the neutral
atom density, unless of course the density matrix is initialized from a previous
calculation by the use of the -i command line parameter; in that case, there
is naturally no point to perform any initial bulk run either, and the parameter
SCF::TwoProbe::InitialBulkRun will be ignored. Experience shows that the
initial bulk run almost always improves the convergence rate and the probability
for successful convergence for two-probe systems.
(Default value: T)
D.2.10 SIMULATION::
These keywords control the simulation steps performed, i.e. geometry relaxation, etc.
Simulation::ConstrainAtoms integer:integer[,integer:integer]
Fixes the positions of certain atoms during a relaxation. Example: a system con-
sists of 12 atoms numbered 0, 1, 2, : : :, 11 (the numbering of atoms always begins
at zero). With the command
Simulation::ConstrainAtoms 0:2,7:10,4:5
we constrain atoms 0, 1, 2, 7, 8, 9, 10, 4, 5. Negative atom indices correspond

to counting backwards from zero, so the command
Simulation::ConstrainAtoms 0:2,-5:-2,-8,5
is equivalent to the one above.

Simulation::Type SingleConfiguration k Relaxation

Parameter defining the overall simulation. A SingleConfiguration simulation
calculates the self-consistent energy and potential, while the ions are moved in
order to minimize the total energy in a Relaxation simulation.
(Default value: SingleConfiguration)
D.2.11 SYSTEM::
Keywords describing the overall symmetry of the system.
System::Type Molecule k Bulk k Electrode k TwoProbe

The symmetry of the system. A Molecule is an isolated system and a Bulk system
is a periodic system. An Electrode system is similar to a bulk system except that
additional information is calculated such that an electrode can be used to define
a two-probe system. Finally, a TwoProbe system an open system where a central
region (the scattering region) is connected with two electrode regions.
D.2.12 TWOPROBE::
This class specifies parameters which are specific to two-probe systems.
TwoProbe::CentralAtomList %block
The coordinates of the atoms in the central region, given in a list (cf. AtomList
D.2 Keywords 101
on page 90).
TwoProbe::LeftElectrode::ATKfile string
Name of the ATK file that defines a calculation of the left electrode. The keyword
can be used instead of the keyword TwoProbe::LeftElectrode::NetCDFfile
in order to avoid a separate calculation of the left electrode.
TwoProbe::LeftElectrode::EquivalentCentralAtom integer
The left electrode is aligned with the central region by requiring that the atom cor-
responding to this parameter is displaced by one lattice vector relative to the atom
corresponding to the parameter TwoProbe::LeftElectrode::EquivalentElectrodeAtom.
If the atom number is negative, it is counted backwards from the end of the atom
list.
(Default value: 0)
TwoProbe::LeftElectrode::EquivalentElectrodeAtom integer
Defines the number of the atom in the left electrode which is aligned with the
atom in the central region specified by TwoProbe::LeftElectrode::EquivalentCentralAtom.
(Default value: 0)
TwoProbe::LeftElectrode::NetCDFfile string
Name of the NetCDF file that defines the left electrode. The keyword can be
substituted by TwoProbe::LeftElectrode::NetCDFfile in order to avoid a sep-
arate calculation of the left electrode.
TwoProbe::LeftElectrode::Voltage real EnergyUnit

Electrochemical potential of the left electrode.
TwoProbe::LeftSurface::NumberOfAtoms integer
Defines which atoms in the central region that are regarded as surface atoms of
the left electrode. These atoms are fixed in relaxations.
(Default value: 0)
TwoProbe::Molecule::GateVoltage real EnergyUnit

Simulate the effect of a gate electrode by shifting the MPSH part of the two-probe
Hamiltonian by this amount. (Default value: 0. eV)
TwoProbe::RightElectrode::ATKfile string
Name of the ATK file that defines a calculation of the right electrode. The
keyword can be used instead of the keyword TwoProbe::RightElectrode::
NetCDFfile in order to avoid a separate calculation of the right electrode.
TwoProbe::RightElectrode::EquivalentCentralAtom integer
Defines the number of the atom in the central region which is used for aligning
the right electrode with the central region. (Default value: -1)
Figure D.1: The contour used for the complex energy integral.
TwoProbe::RightElectrode::EquivalentElectrodeAtom integer
Defines the number of the atom in the right electrode which is aligned with the
atom in the central region specified by TwoProbe::RightElectrode::EquivalentCentralAtom.
(Default value: -1)
TwoProbe::RightElectrode::NetCDFfile string
Name of the NetCDF file which defines the right electrode. The keyword can
be substituted by TwoProbe::RightElectrode::NetCDFfile in order to avoid a
separate calculation of the right electrode.
TwoProbe::RightElectrode::Voltage real EnergyUnit

Electrochemical potential of the right electrode.
TwoProbe::RightSurface::NumberOfAtoms integer
Defines which atoms in the central region that are regarded as surface atoms of
the right electrode. These atoms are fixed in relaxations.
(Default value: 0)
TwoProbe::UseMultigridForElectrostatics logical
When set to T, the Hartree potential is solved with a multigrid method, and this
allows for different left and right electrodes. In this case it is crucial that the size
of the unit cells in the A and B directions are the same for both electrodes (i.e.
not the same for A and B, but the same in the A direction in the left and right
electrodes, and the same for the B direction). When this option is set to F, the
Hartree potential is solved using FFT and periodic boundary conditions, and in
this case the left and right electrodes must be the same.
(Default value: F)
TWOPROBE::CONTOUR
This subclass of keywords is related to the contour integration.

D.2 Keywords 103
TwoProbe::Contour::Circle::NumPoints integer
The number of points along the arc part of the complex contour starting at
TwoProbe::Contour::DeltaEMin and ending at the Fermi energy Ef .
(Default value: 30)
TwoProbe::Contour::DeltaEMin real EnergyUnit

The distance from the Fermi energy Ef to the starting energy of the complex energy
contour. In a ATK run this value should be larger than the binding energy of the
lowest eigenstate, otherwise there will be poles outside the contour and charge
will be missing in the integration.
TwoProbe::Contour::FermiLine::NumPoints integer
This keyword determines the number of points from the Fermi energy Ef to 1. The
Fermi function is used as a weight function, and the important energy range has a
width of about 4kB T where kB is the Boltzmann constant and T the temperature.
(Default value: 10)
TwoProbe::Contour::KTPoles::NumPoints integer
The complex contour always stays a distance from the real axis. This means
that it encloses some poles of the Fermi functions. In ATK the number of poles
is specified by this keyword, and this indirectly determines the distance of the
contour from the real axis.
(Default value: 4)
TwoProbe::Contour::RealAxis::PointDensity real EnergyUnit

The spacing between points on the close-to-real axis part of the contour in the
voltage bias window.
(Default value: 0.02 eV)
TwoProbe::Contour::RealAxis::Eta real EnergyUnit

Small complex number added to the denominator in the real axis integration.
(Default value: 0.01 eV)
D.2.13 UNITCELL::
This class defines the unit cell of the calculation.
UnitCell::LatticeConstant real LengthUnit

The lattice constant of the unit cell.
(Default value: 1. Ang)
UnitCell::X0 real LengthUnit

Cartesian x coordinate of the starting point for the cell vectors.
UnitCell::Y0 real LengthUnit

Cartesian y coordinate of the starting point for the cell vectors.
UnitCell::Z0 real LengthUnit

Cartesian z coordinate of the starting point for the cell vectors.
A P P E N D I X
E
OUTPUT FILE
During a run with ATK, a log or output file is always generated. In addition to informing
the user of the progress of the calculations, it also contains many important calculated
quantities such as total energy or final positions of a relaxation run. The log file also
documents the parameters used in the calculation.
Most parts of the log file are self-documented and easy enough to understand just by
looking at it. Here we will just point out some details which perhaps are not so obvious,
usually due to special notation.
E.1 HEADER AND SELF-CONSISTENT CYCLE
The first portion of the log file, after the header which tells when, where, how and by
whom ATK was launched, is simply a record of the parameters used in the calculations.
Parameters not found in the input atk file are taken from the defaults file, but the log file
does not become a copy of the defaults file, since only parameters which are actually
relevant for the particular kind of calculation desired are displayed.
Unless ATK was started in analysis mode (using the -a option), the next part of the log
file is almost always the progress information from the self-consistent cycle. Except for
the first step, this consists of a single line of the form (it is all a single line in the log file)
# sc 1 : FermiEnergy = -3.8255 eV Ebs = -126.9 eV

dRho = 2.2543E-2 dEbs = 1.43532E-4 dH = 6.34533E-01
Here, "Ebs" is the band structure energy, i.e. the sum of the energies of all occupied
states, while "dRho", "dEbs" and "dH" are the relative changes in the density matrix,
band structure energy, and Hamiltonian, respectively, compared to the previous step.
105
106 Output file
These are the quantities which are compared to the tolerance, set by SCF::Tolerance,
which must be met for convergence.
For a two-probe run with open boundary conditions (there is also, if so requested, an
initial run with periodic boundary conditions, as explained above in relation to the
keyword SCF::TwoProbe::InitialBulkRun), the line is slightly different:
# sc 3 : q = 8.13734 e Ebs = -20.45441 eV

dRho = 4.0454E-02 dEbs = 3.3454E-03 dH = 4.3454E-02
Here "q" is the total charge in the central region, in units of the electron charge. For
molecular, bulk and electrode systems, the total charge is of course strictly conserved,
but this is not the case in a two-probe system (in particular under bias), where we have
open boundary conditions.
Note carefully that the band structure energy "Ebs" is defined in a different way in the
initial bulk run compared to a two-probe system with open boundary conditions, which
is why this quantity will exhibit a jump when going from one self-consistent cycle to
the other. In the first case, the electrodes are included, in the latter they are not.
When running with the command line option --verbose, there is additional progress
information for each step in the self-consistent cycle. First of all, the Mulliken popula-
tions are displayed for each step (this can often be a useful piece of information when
trying to cure convergence problems), and in addition there is a total energy summary
in each step.
Moreover, with the --verbose option, ATK also prints the k-point sampling grid used
for the calculation to the output file. The output looks something like
%block BrillouinZoneSampling
# ----------------------------------------------------------------
# KX (1/Ang) KY (1/Ang) KZ (1/Ang) Weight
# ----------------------------------------------------------------
-0.38500 -0.38500 -0.38500 0.25000
-1.15500 0.38500 0.38500 0.25000
0.38500 -1.15500 0.38500 0.25000
-0.38500 -0.38500 1.15500 0.25000
# ----------------------------------------------------------------
%endblock BrillouinZoneSampling
for a 222 grid in fcc Au with lattice constant 4.08 Å. Note that 222=8, but 4 of
these can be related to the other 4 by time-reversal symmetry, which is used by ATK to
reduce the number of required calculations.
E.2 Total energy and atom list 107
E.2 TOTAL ENERGY AND ATOM LIST
The total energy information is always displayed after the self-consistent loop has com-
pleted. It is a large block,
%block Results::TotalEnergy
.
.
.
%endblock Results::TotalEnergy
which consists of two main parts. First, there is a summary of all the different terms in
the total energy:
# Total Energy = -465.33332 eV

# Potential Energy = -465.33332 eV
# Ion Kinetic Energy = 0.00000 eV
# --------------------------------------
# NumOrb: Energy of system
# --------------------------------------
# Charge = 8.0000
# FermiEnergy = -0.7348 eV
# Band Structure Energy = -101.6476 eV
# --------------------------------------
# El. Kinetic Energy = 337.3125 eV
# NL Energy = -113.9817 eV
# dEna = -5.4525 eV
# dUscf = 1.0067 eV
# XC Energy = -116.0787 eV
# ENA = 69.3998 eV
# Eions = -637.5394 eV
# --------------------------------------
# Electron kinetic energy = 337.3125 eV
# Electrostatic energy = -686.5671 eV
# Exchange-corr. energy = -116.0787 eV
# --------------------------------------
# Total potential energy = -465.3333 eV
# --------------------------------------
The first line, the total energy, is one of the most important parameters in a DFT cal-
culation. The total energy (of the ions, as it is) is the sum of the ionic potential energy
and the kinetic energy. For the calculations that presently can be performed in ATK the
ions are stationary, and thus the total and potential energies are always equal.
108 Output file
Note carefully that since ATK uses pseudopotentials to represent the core electrons,
the absolute energies reported in the output files have no real meaning. Only energy
differences have physical relevance, be they between different states in the same sys-
tem, or the same system calculated for different conditions (such as when relaxing the
geometry).
The total energy block also reports the charge in the system, and the Fermi energy and
band structure energy. As mentioned already, the band structure energy is the sum of
the energies of all eigenstates, weighted with their occupation, and the total charge in
the system is conserved for molecular, bulk and electrode systems. Note that for two-
probe systems, the charge listed in the total energy block includes the charges located
in the two electrode unit cells as well, whereas the value of "q" as presented during the
self-consistent cycle only counts the charges in the central region.
For bulk and molecule calculations, the Fermi energy is calculated by requiring the
system to be charge neutral. More precisely, the Fermi level is the energy which makes
the Fermi occupation of the orbitals equal to the number of electrons in the system.
Note that for molecules, if there is a HOMO–LUMO gap, the Fermi level can vary
dramatically as a function of the temperature and other parameters, and is thus not
really well-defined. It is then physically more interesting to look at the HOMO level.
The Fermi energy is not presented for two-probe systems, since it is an equilibrium
quantity which is not well-defined when open boundary conditions apply.
The total energy is made up of contributions from several terms, which are also pre-
sented in this block. For a more detailed presentation of the definitions of these, please
refer to Ref. 11. Finally, the total potential energy (or more specifically the potential
energy of the ions) is the sum of the ground state energy of the electronic system and
the ion–ion interaction, and can be subdivided into three terms: the electron kinetic
energy (i.e. the kinetic energy of the non-interacting electron gas), the electrostatic
energy (the mean-field electrostatic energy of the electron gas and the ions), and the
exchange–correlation energy of the electron gas.
The second part of the total energy block is a report on the atomic configuration, listing
the positions (i.e. the final positions in the case of a relaxation run), plus the velocities
and forces on each atom:
# ----------------------------------------------------------------
# Number Of Atoms = 3
# ----------------------------------------------------------------
# Positions in Angstrom
# Velocities in Angstrom/fs
# Forces in eV/Ang
# ----------------------------------------------------------------
# X Y Z FX FY FZ
# ----------------------------------------------------------------
0.757 0.586 0.000 5.419E-01 5.080E-01 7.978E-14
-0.757 0.586 0.000 -5.419E-01 5.080E-01 1.153E-13
E.3 Mulliken population 109
0.000 0.000 0.000 -5.608E-14 -1.017E+00 5.326E-14

# ----------------------------------------------------------------
# Total Force (eV/Ang) = 0.0 -0.0014080 0.0
# Total Magnitude (eV/Ang) = 0.00140804
# Largest force component (eV/Ang) = 1.01733
# Root-mean-squared Force Magnitude (eV/Ang) = 0.844244
# ----------------------------------------------------------------
%endblock Results::TotalEnergy
Actually, the example above is a slightly modified version of a real output from ATK. We
have removed the velocities, which naturally always zero are in the presently supported
calculations, simply in order to fit the output.
Unless you are running a relaxation, it is recommended to always check in the output
that the forces are not too large. If they are, it is an indication that the system setup
does not correspond to a natural configuration, and the results may be questionable for
this reason. The forces calculated by ATK are analytical DFT forces [11].
During a relaxation, there will be one self-consistent cycle per iteration in the conjugent
gradient minimization of the total energy. Each step will display the self-consistent total
energy block, plus the atomic configuration in that step, so that one can follow the
progress of how the atoms are moved into their relaxed positions.
E.3 MULLIKEN POPULATION
If the keyword Analysis::PrintMullikenPopulation is set to true (this is the default),

and if ATK is run with the command line option --verbose, the Mulliken populations
will be evaluated and printed to the log file for every step in the self-consistent cycle, as
well as at the end of the calculation. Unless one is having trouble with the convergence
of the self-consistent cycle, however, it is often not necessary to use the --verbose
option, since it bloats the output file quite a bit. In that case, it is possible to perform
the Mulliken population analysis as a simple post-processing analysis step, by running
atk file.atk mulliken.out -a file.nc --verbose --force
where file.atk is the name of the original input file for the calculation (which natu-
rally must not set the keyword Analysis::PrintMullikenPopulation to False!), and
file.nc is the NETCDF file produced by the original run. The --force option is re-
quired to overwrite file.nc, but in fact the file will not be modified.
The Mulliken populations are presented in a block with, except for a report on the total
charge, simply consists of one line per atom in the setup (for two-probes, for each atom
in the central region). For each atom, the total Mulliken population, which is the same
110 Output file
` m Label Solid harmonic ` m Label Solid harmonic
0 0 s 1 3 –3 f(–3) y(3x2 y2 )
3 –2 f(–2) xyz
1 –1 py y 3 –1 f(–1) y(4z 2 x2 y2 )
1 0 pz z 3 0 f(0) z (2z 2 3x2 3y2 )
1 +1 px x 3 +1 f(1) x(4z 2 x2 y2 )
3 +2 f(2) z (x2 y2 )
2 –2 d(–2) xy 3 +3 f(3) x(x2 3y 2 )
2 –1 d(–1) yz
2 0 d(0) 2z 2 (x2 + y2 )
2 +1 d(1) xz
2 +2 d(2) x2 y2
Table E.1: ATK employs solid harmonics to expand the angular part of the basis orbitals. Note
that normalization factors are omitted from all expressions here.
as the total number of electrons on that atom, is listed as "q", and after that follows
a breakdown of how these electrons are distributed among the basis orbitals. These
orbitals are labelled as follows
An example of such output, where a DZP basis set was used for both H and Au, could
look like this:
%block Results::MullikenPopulations
0 H q = 1.00000 [ s: 1.014, s: -0.014, y: 0.000, z: 0.000, x: 0.000 ]
1 Au q = 11.00000 [ s: -0.006, s: 2.006, d(-2): 1.840, d(-1): 1.840,
d(0): 1.331, d(1): 1.840, d(2): 1.331, d(-2): 0.160, d(-1): 0.160,
d(0): 0.169, d(1): 0.160, d(2): 0.169, y: 0.000, z: 0.000, x: 0.000 ]
# ----------------------------------------------------------------
# Total charge = 12
# ----------------------------------------------------------------
%endblock Results::MullikenPopulations
There are two "s" orbitals for H, which is the "DZ" part of the basis set, plus the "P"
(polarization) "p" orbitals. For Au, we have two sets of "s" and "d" orbitals (again "DZ"),
plus a single set of polarization "p" orbitals. Note that the valence configuration for Au
is 5d10 6s2 , and thus the next higher shell, i.e. the polarization shell, is 6p.
E.4 Spectra 111
E.4 SPECTRA
E.4.1 ENERGY SPECTRUM
For molecular systems, the last quantity reported in the output file of a regular calcula-
tion is the molecular eigenspectrum (provided the keyword Analysis::PrintMolecularEigenvalues
is set to true, which is the default value). The energy zero-level is chosen at the vacuum
level, i.e. the value of the effective potential far away from the molecule.
In addition, the energies of the HOMO and LUMO levels are presented, and the oc-
cupation of each eigenstate is indicated. Note that the HOMO and LUMO are not
entirely well-defined for spin-dependent systems.
For bulk systems, the vacuum level is not well defined, in this case the energy zero
is arbitrary. Instead, the energy spectrum (i.e. the output obtained from specifying
Analysis::BandLine keywords) is always presented relative the Fermi energy. All
wave vectors in the output file are reported in inverse Ångström.
Finally, for two-probe systems, the MPSH eigenstates are reported relative to the aver-
age Fermi level of the electrodes. Thus, if an equal voltage (say, 1 V) is applied to both
electrodes, the MPSH spectrum does not change.
E.4.2 TRANSMISSION SPECTRUM
The last part of a two-probe output file contains the energy-dependent transmission
spectrum, the details of which are controlled by the keywords Analysis::TransmissionSpectrum.
As for the MPSH eigenstates, the energy scale for the transmission spectrum is reported
relative the average Fermi level.
Note that the transmission coefficient can be larger than unity if multiple transmission
channels are available at a certain energy. Also, when the transmission coefficient is
less than unity, it does not mean that this state is not transmitted; it only means that
the probability of the incoming wave function to be transmitted is smaller than one
(quantum mechanics is inherently statistical).
If the keywords Analysis::TransmissionSpectrum::NumKPoints::A/B are different

from unity (and in general they should be, in order to obtain a converged transmission
spectrum; the most notable exception being that of effectively one-dimensional elec-
trodes), the transmission is first given as a weighted average (i.e. the actual transmission
spectrum), and then the spectrum for each sampled k-point in the 2D surface Brillouin
zone is reported.
The two last column of the (average) transmission spectrum are the total density of
112 Output file
states (DOS) of the entire two-probe system, and the partial density of states (PDOS)
of the central region excluding surface atoms, defined by the parameters TwoProbe::
LeftSurface::NumberOfAtoms and TwoProbe::RightSurface::NumberOfAtoms.
For the k-point resolved spectra, there are also additional columns giving the largest
eigenvalues of the decomposed transmission eigenchannels, as controlled by Analysis::
TransmissionSpectrum::NumEigenChannels.
Finally, the output file also contains the calculated current (at finite bias) and conduc-
tance. The value reported in the output file is essentially the transmission at the Fermi
level (times the quantum conductance unit, 77.48 S); to calculate the differential con-
ductance, it is necessary to perform a numerical differentiation of the current–voltage
function.
A P P E N D I X
F
WRITING SCRIPTS
Because ATK takes a number of command line options, writing small scripts is a very
powerful and yet simple way to perform complex calculations, as described in this
tutorial.
One may create these scripts in any suitable scripting language which is able to call
other executables. In the following sections we provide examples of such scripts
for Linux/Cygwin (bash shell scripts) and Windows (DOS batch files), plus platform-
independent scripts written in Python. All examples are related to tasks considered in
the quick tour section of this tutorial, as referred to in each case below.
As with the input files it is recommended that you type and run the scripts yourself, but
for reference they can also be found in the examples/tutorial/quicktour directory.
F.1 BASH SCRIPTS FOR LINUX AND CYGWIN
In order to execute the scripts you write, it is necessary to make them executable
through the command
$ chmod +x script.sh
where script.sh of course is the name of the script and the $ sign symbolizes the
command prompt. To run the script, simply type
$ ./script.sh
at the command prompt.
113
114 Writing Scripts
Note! In Linux (and similar environments such as Cygwin) the

command interpreter or shell interprets each line in a script file as
a separate command. However, many of the commands we wish
to run are rather long, and do not fit on a single line. For improved
readability we have therefore in the script listings below used a "n"
sign at the end of the line. This is a "continuation sign", which
tells the command shell that the text on the following line is to be
considered as a continuation of the present one, and not as a new
command. If you wish, you may of course omit the continuation
sign and enter the entire command on a single line.
The following script (latenergy.sh) calculates the total energy for a series of lattice
constants for the lithium chain discussed on page 10:
#!/bin/sh
for lat in 2.82 2.84 2.86 2.88 2.90 2.92 ; do

atk lichain.atk lichain$lat.out \
--force -o lichain$lat.nc \
--param "UnitCell::LatticeConstant $lat Angstrom" $*
done
Next, we have a script ivrun.sh that calculates the current–voltage characteristics of

the Li–H–H–Li two-probe system considered on page 17:
#!/bin/sh
oldbias=0.0
for bias in 0.2 0.4 0.6 0.8 1.0 ; do

atk lih2-0.0.atk lih2-$bias.out \
--force -o lih2-$bias.nc -i lih2-$oldbias.nc \
--param "TwoProbe::LeftElectrode::Voltage $bias eV" $*
oldbias=$bias
done
For the same system, we also made an I–V calculation when the geometry was relaxed
on page 19, using the script ivrunrelax.sh:
#!/bin/sh
oldbias=0.0
bias=0.0
atk lih2-$bias-relax.atk lih2-$bias-relax.out \

-o lih2-$bias-relax.nc --force $*
F.2 Batch Files for Windows 115
for bias in 0.0 0.2 0.4 0.6 0.8 1.0 ; do

atk lih2-0.0-relax.atk lih2-$bias-relax.out \
-o lih2-$bias-relax.nc -i lih2-$oldbias-relax.nc \
--param "TwoProbe::LeftElectrode::Voltage $bias eV" \
--force $*
oldbias=$bias
done
Note how the calculation for the previous value of the bias is used as an initial guess
for the density matrix by using the command -i.
Again for the Li–H–H–Li, we also considered the influence of a gate electrode (page 20).
The corresponding script gaterun.sh is:
#!/bin/sh
oldgate=0.0
cp lih2-0.0.nc gate-0.0.nc
for gate in 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 ; do
atk lih2-0.0.atk gate-$gate.out \
-o gate-$gate.nc -i gate-$oldgate.nc \
--param "TwoProbe::Molecule::GateVoltage -$gate eV" \
--force $*
oldgate=$gate
done
Note that we use a negative gate voltage, but in order to avoid strange filenames we in-
clude the negative sign explicitly in the --param command and loop over the absolute
value of the potential.
F.2 BATCH FILES FOR WINDOWS
In Windows, script files are usually referred to as batch files, and are given the exten-
sion .bat. Such files can be executed directly by double-clicking them in Windows
Explorer.
116 Writing Scripts
Note! In Windows it is not possible to break long lines in a batch

file; the commands must be given on a single line. However, many
of the commands we wish to run are rather long, and do not fit on a
single line. For improved readability we have therefore in the script
listings below anyway broken them, by using a "n" sign at the end
of the line. This character must be removed, and the corresponding
lines joined together, if you wish to copy/paste the batch files from
this manual. The batch files provided in the examples directory are
however correctly formatted already.
The following script (latenergy.bat) calculates the total energy for a series of lattice
constants for the lithium chain discussed on page 10 (this must all be entered on a
single line, either in a batch file or on a Windows command line; c.f. the box on
page 116):
for %%A in (2.82 2.84 2.86 2.88 2.90 2.92) do \

atk lichain.atk lichain%%A.out -o lichain%%A.nc \
--param "UnitCell::LatticeConstant %%A Ang" --force
Next, we have a script ivrun.bat that calculates the current–voltage characteristics of

set oldbias=0.0
for %%V in (0.2 0.4 0.6 0.8 1.0) do call :Loop %%V
goto eof
:Loop
atk lih2-0.0.atk lih2-%1.out -o lih2-%1.nc -i lih2-%oldbias%.nc \
--param "TwoProbe::LeftElectrode::Voltage %1 eV" --force
set oldbias=%1
:eof
on page 19, using the script ivrunrelax.bat:
set oldbias=0.0
atk lih2-0.0-relax.atk lih2-0.0-relax.out -o lih2-0.0-relax.nc --force
for %%V in (0.2 0.4 0.6 0.8 1.0) do call :Loop %%V
goto eof
F.3 Python Scripts 117
:Loop
atk lih2-0.0-relax.atk lih2-%1-relax.out -o lih2-%1-relax.nc \
-i lih2-%oldbias%-relax.nc --force \
--param "TwoProbe::LeftElectrode::Voltage %1 eV"
set oldbias=%1
:eof
The corresponding script gaterun.bat is:
set oldgate=0.0
copy /Y lih2-0.0.nc gate-0.0.nc
for %%A in (1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0) do call :Loop %%A
goto eof
:Loop
atk lih2-0.0.atk gate-%1.out -o gate-%1.nc -i gate-%oldgate%.nc \
--param "TwoProbe::Molecule::GateVoltage -%1 eV" --force
oldgate=%1
:eof
F.3 PYTHON SCRIPTS
In order to use Python scripts, you must have the Python interpreter installed on your
computer. To test if you do, simply type the command
$ python -V
at the command prompt in either Windows or Linux. If Python is installed, this com-
mand should display the Python version number. Be careful to type a capital V; if you
type python -v you will enter the interactive mode in Python, and you must press
Ctrl+D to exit.
Python scripts are executed by simply giving the script name as the only argument to
Python:
$ python script.py
118 Writing Scripts
where script.py is the name of the script.
The following script (latenergy.py) calculates the total energy for a series of lattice
constants for the lithium chain discussed on page 10:
import os
latconst = ["2.82", "2.84", "2.86", "2.88", "2.90", "2.92"]
for lat in latconst:

os.popen("atk lichain.atk lichain"+lat+".out \
--param \"UnitCell::LatticeConstant "+lat+" Ang\" \
-o lichain"+lat+".nc" --force)
Next, we have a script ivrun.py that calculates the current–voltage characteristics of

import os
bias = ["0.2", "0.4", "0.6", "0.8", "1.0"]

oldV = 0.0
for V in bias:
os.popen("atk lih2-0.0.atk lih2-"+V".out \
--param "TwoProbe::LeftElectrode::Voltage "+V" eV" \
--force -i lih2-"+oldV".nc -o lih2-"+V+".nc")
oldV = V
on page 19, using the script ivrunrelax.py:
import os
bias = ["0.2", "0.4", "0.6", "0.8", "1.0"]

oldV = 0.0
V = 0.0
os.popen("atk lih2-"+V"-relax.atk lih2-"+V"-relax.out \

-o lih2-"+V"-relax.nc" --force)
for V in bias:
os.popen("atk lih2-0.0-relax.atk lih2-"+V"-relax.out \
--param "TwoProbe::LeftElectrode::Voltage "+V" eV" \
-i lih2-"+oldV"-relax.nc -o lih2-"+V+"-relax.nc" --force)
oldV = V
F.3 Python Scripts 119
The corresponding script gaterun.py is:
import os
os.popen("cp lih2-0.0.nc gate-0.0.nc")
gatepotential = ["1.0", "2.0", "3.0", "4.0", "5.0", \

"6.0", "7.0", "8.0", "9.0", "10.0"]
oldVg = 0.0
for Vg in gatepotential:
os.popen("atk lih2-0.0.atk gate-"+Vg".out \
--param "TwoProbe::Molecule::GateVoltage -"+Vg" eV" \
--force -i gate-"+oldVg".nc -o gate-"+Vg+".nc")
oldVg = Vg
120 Writing Scripts
A P P E N D I X
G
ATOMIC DATA
This appendix list the pseudopotentials included in ATK 2.0, along with the corre-
sponding atomic configurations used for constructing the basis sets. From the valence
charge of each pseudopotential, the valence configuration can be obtained. In the Sin-
gle Zeta (SZ) basis set there is one basis function for each valence orbital. For instance,
for a Au (Z=79) pseudopotential with ionic charge 11, the valence orbitals will consist
of 5d10 6s1 .
The polarization orbitals will have lowest angular momentum not represented in the
SZ basis set. For instance, for the Au system, the polarization orbital will be a p orbital.
Z Symbol Electron configuration Valence configuration Valence charge
1 H 1s1 1s1 1
2 He 1s 2 2s2 2
3 Li [He]2s1 2s1 1
4 Be [He]2s 2 2s2 2
5 B [He]2s2 2p1 2s2 2p1 3
6 C 2
[He]2s 2p 2 2s2 2p2 4
7 N [He]2s2 2p3 2s2 2p3 5
8 O 2
[He]2s 2p 4 2s2 2p4 6
9 F 2
[He]2s 2p 5 2s2 2p5 7
10 Ne [He]2s2 2p6 2s2 2p6 8
11 Na [Ne]3s 1 3s1 1
12 Mg [Ne]3s2 3s2 2
continued on next page
121
122 Atomic Data
13 Al [Ne]3s2 3p1 3s2 3p1 3

14 Si [Ne]3s2 3p2 3s2 3p2 4
15 P 2
[Ne]3s 3p 3 3s2 3p3 5
16 S [Ne]3s2 3p4 3s2 3p4 6
17 Cl 2
[Ne]3s 3p 5 3s2 3p5 7
18 Ar [Ne]3s2 3p6 3s2 3p6 8
19 K 0
[Ar]3d 4s 1 4s1 1
20 Ca [Ar]3d0 4s2 4s2 2
21 Sc 1
[Ar]3d 4s 2 4s2 3d1 3
22 Ti 2
[Ar]3d 4s 2 4s2 3d2 4
23 V [Ar]3d3 4s2 4s2 3d3 5
24 Cr 4
[Ar]3d 4s 2 4s2 3d4 6
25 Mn [Ar]3d5 4s2 4s2 3d5 7
26 Fe 6
[Ar]3d 4s 2 4s2 3d6 8
27 Co [Ar]3d7 4s2 4s2 3d7 9
28 Ni 8
[Ar]3d 4s 2 4s2 3d8 10
29 Cu 10
[Ar]3d 4s 1 4s1 3d10 11
30 Zn [Ar]3d10 4s2 4s2 3d10 12
31 Ga 10
[Ar]3d 4s 4p2 1 4s2 4p1 3
32 Ge [Ar]3d10 4s2 4p2 4s2 4p2 4
33 As 10
[Ar]3d 4s 4p2 3 4s2 4p3 5
34 Se [Ar]3d10 4s2 4p4 4s2 4p4 6
35 Br 10
[Ar]3d 4s 4p2 5 4s2 4p5 7
36 Kr [Ar]3d10 4s2 4p6 4s2 4p6 8
37 Rb 0
[Kr]4d 5s 1 5s1 1
38 Sr 0
[Kr]4d 5s 2 5s2 2
39 Y [Kr]4d1 5s2 5s2 4d1 3
40 Zr 2
[Kr]4d 5s 2 5s2 4d2 4
41 Nb [Kr]4d4 5s1 5s2 4d3 5
42 Mo 5
[Kr]4d 5s 1 5s1 4d5 6
43 Tc [Kr]4d5 5s2 5s2 4d5 7
44 Ru 7
[Kr]4d 5s 1 5s2 4d6 8
45 Rh 8
[Kr]4d 5s 1 5s1 4d8 9
46 Pd [Kr]4d10 5s0 4d10 10
47 Ag 10
[Kr]4d 5s 1 5s1 4d10 11
48 Cd [Kr]4d10 5s2 5s2 4d10 12
49 In 10 2
[Kr]4d 5s 5p 1 5s2 5p1 3
50 Sn [Kr]4d10 5s2 5p2 5s2 5p2 4
51 Sb 10 2
[Kr]4d 5s 5p 3 5s2 5p3 5
123
52 Te [Kr]4d10 5s2 5p4 5s2 5p4 6

53 I [Kr]4d10 5s2 5p5 5s2 5p5 7
54 Xe 10 2
[Kr]4d 5s 5p 6 5s2 5p6 8
55 Cs [Xe]5d0 6s1 6s1 5p6 7
56 Ba 0
[Xe]5d 6s 2 6s2 5p6 8
57 La [Xe]5d1 6s2 6s2 5d1 3
58 Ce 2 0
[Xe]4f 5d 6s 2 6s2 4f2 4
59 Pr [Xe]4f3 5d0 6s2 6s2 4f3 5
60 Nd 4 0
[Xe]4f 5d 6s 2 6s2 4f4 6
61 Pm 5 0
[Xe]4f 5d 6s 2 6s2 4f5 7
62 Sm [Xe]4f6 5d0 6s2 6s2 4f6 8
63 Eu 7 0
[Xe]4f 5d 6s 2 6s2 4f7 9
64 Gd [Xe]4f7 5d1 6s2 6s2 5d1 4f7 10
65 Tb 9 0
[Xe]4f 5d 6s 2 6s2 4f9 11
66 Dy [Xe]4f10 5d0 6s2 6s2 4f10 12
67 Ho 11
[Xe]4f 5d 6s0 2 6s2 4f11 13
68 Er 12
[Xe]4f 5d 6s0 2 6s2 4f12 14
69 Tm [Xe]4f13 5d0 6s2 6s2 4f13 15
70 Yb 14
[Xe]4f 5d 6s0 2 6s2 4f14 16
71 Lu [Xe]4f14 5d1 6s2 6s2 5d1 4f14 17
72 Hf 14
[Xe]4f 5d 6s2 2 6s2 5d2 4
73 Ta [Xe]4f14 5d3 6s2 6s2 5d3 5
74 W 14
[Xe]4f 5d 6s4 2 6s1 5d5 6
75 Re [Xe]4f14 5d5 6s2 6s2 5d5 7
76 Os 14
[Xe]4f 5d 6s6 2 6s2 5d6 8
77 Ir 14
[Xe]4f 5d 6s7 2 6s2 5d7 9
78 Pt [Xe]4f14 5d10 6s0 6s1 5d9 10
79 Au 14
[Xe]4f 5d 6s10 1 6s1 5d10 11
80 Hg [Xe]4f14 5d10 6s2 6s2 5d10 12
81 Tl 14 10
[Xe]4f 5d 6s 6p 2 1 6s2 6p1 3
82 Pb [Xe]4f14 5d10 6s2 6p2 6s2 6p2 4
83 Bi 14 10
[Xe]4f 5d 6s 6p 2 3 6s2 6p3 5
84 Po 14 10
[Xe]4f 5d 6s 6p 2 4 6s2 6p4 6
85 At [Xe]4f14 5d10 6s2 6p5 6s2 6p5 7
86 Rn 14 10
[Xe]4f 5d 6s 6p 2 6 6s2 6p6 8
87 Fr [Rn]7s1 7s1 6p6 7
88 Ra [Rn]7s 2 7s2 6p6 8
89 Ac [Rn]6d1 7s2 5f0 7s2 6d1 3
90 Th 2 2
[Rn]6d 7s 5f 0 7s2 6d2 4
124 Atomic Data
91 Pa [Rn]6d1 7s2 5f2 7s2 6d1 5f2 5

92 U [Rn]6d1 7s2 5f3 7s2 6d1 5f3 6
93 Np [Rn]6d0 7s2 5f5 7s2 6d1 5f4 7
94 Pu [Rn]6d0 7s2 5f6 7s2 5f6 8
95 Am [Rn]6d0 7s2 5f7 7s2 5f7 9
96 Cm [Rn]6d1 7s2 5f7 6d1 7s2 5f7 10
97 Bk [Rn]6d0 7s2 5f9 7s2 5f9 11
98 Cf [Rn]6d0 7s2 5f10 7s2 5f10 12
99 Es [Rn]6d0 7s2 5f11 7s2 5f11 13
100 Fm [Rn]6d0 7s2 5f12 7s2 5f12 14
101 Md [Rn]6d0 7s2 5f13 7s2 5f13 15
102 No [Rn]6d0 7s2 5f14 7s2 5f14 16
103 Lr [Rn]6d1 7s2 5f14 7s2 6d1 5f14 17
References
[1] M. BRANDBYGE, J.-L. MOZOS, P. ORDEJÓN, J. TAYLOR, AND K. STOKBRO, Density-

functional method for nonequilibrium electron transport, Phys. Rev. B 65, 165401
(2002). i, 46
[2] C. G. BROYDEN, A class of methods for solving nonlinear simultaneous equations,

Math. Comput. 19, 577–593 (1965). 24
[3] W. H. BUTLER, X.-G. ZHANG, T. C. SCHULTHESS, AND J. M. MACLAREN, Spin-

dependent tunneling conductance of fe|mgo|fe sandwiches, Phys. Rev. B 63,
54416 (2001). 59
[4] D. M. CEPERLEY AND B. J. ALDER, Ground state of the electron gas by a stochastic
method, Phys. Rev. Lett. 45, 566 (1980). 27, 96
[5] P. HOHENBERG AND W. KOHN, Inhomogeneous electron gas, Phys. Rev. 136, 864
(1964). 2
[6] H. J. MONKHORST AND J. D. PACK, Special points for Brillouin-zone integrations,

Phys. Rev. B 13, 5188 (1976). 36, 91, 92
[7] J. P. PERDEW, K. BURKE, AND M. ERNZERHOF, Generalized gradient approximation

made simple, Phys. Rev. Lett. 77, 3865 (1996). 2, 27, 96
[8] J. P. PERDEW AND A. ZUNGER, Self-interaction correction to density-functional

approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981). 2, 23,
27, 96
[9] P. PULAY, Chem. Phys. Lett. 73, 393 (1980). 24
[10] R. H. M. SMIT, Y. NOAT, C. UNTIEDT, N. D. LANG, M. C. VAN HEMERT, AND

J. M. VAN RUITENBEEK, Measurement of the conductance of a hydrogen molecule,
Nature 419, 906 (2002). 5, 17
[11] J. M. SOLER, E. ARTACHO, J. D. GALE, A. GARCÍA, J. JUNQUERA, P. ORDEJÓN, AND

D. SÁNCHEZ-PORTAL, The SIESTA method for ab initio order-N materials simula-
tion, J. Phys. Condens. Matter 14, 2745 (2002). i, 75, 108, 109
[12] K. STOKBRO, J. L. MOZOS, P. ORDEJÓN, M. BRANDBYGE, AND J. TAYLOR, Theoret-

ical study of the nonlinear conductance of Di-thiol benzene coupled to Au(111)
surfaces via thiol and thiolate bonds, Comp. Mat. Sci. 27, 151 (2003). 46
125
[13] J. TAYLOR, H. GUO, AND J. WANG, Ab inito modeling of quantum transport prop-
erties of molecular electronic devices, Phys. Rev. B 63, 245407 (2001). i
[14] Y. TSERKOVNYAK, A. BRATAAS, G.E.W. BAUER, AND B. I. HALPERIN, Cond-

mat/0409242v2 (2005). 59
[15] D. WORTMANN, G. BIHLMAYER, AND S. BLÜGEL, ab initio calculations of interface

effects in tunneling through MgO barriers on Fe(100), J. Phys. Cond. Mat. 16,
S5819 (2004). 60
Index
.nc file, see NetCDF file EnergyWindow, 89

.out file, see output file MPSH::X::, 17
# (comment symbol), 36, 84 Index, 89
$ (command prompt), 3 Label, 89
Spin, 89
A/B/C directions, 8, 36, 43 PrintMolecularEigenvalues, 89, 111
accuracy, 2, 27–28 PrintMullikenPopulation, 89, 109
basis set, 28 TransmissionSpectrum::, 15, 45, 49
force tolerance, 26 Calculate, 89
analysis, 7, 15, 17, 19 E0, 90
Analysis::, 7, 15 E1, 90
3DGrids:: Eta, 47, 90
SaveEffectivePotential, 71 NumEigenChannels, 90, 112
SaveElectronDensity, 71 NumKPoints, 90
SaveTo, 71 NumPoints, 90
BandLine::NumPoints, 87 approximations, 2
BandLine::X::, 73 ATK_DATA_DIR, 81
End, 37, 87 AtomList::
NumPoints, 87 Format, 6, 36, 42, 91
Start, 37, 87 Scale, 36, 91
BlochState::X::
band structure, 11, 36
Index, 87
energy, 105
KPoint, 87
basis set, 9, 23, 28
Label, 87
bias, 17, 19, 50
Spin, 87
boundary conditions, 11, 15, 29, 43
DOS::X::, 17
Broyden mixing, 24
Energy, 87
bulk systems, 7–11
Kpoint, 88
Bulk::
Label, 88
AtomList, 36, 91
Spin, 88
ElectronTemperature, 39, 91
Type, 88
FermiTemperature, 72
MolecularOrbital::X::, 7
NumKPoints, 8, 36, 39, 91
Index, 88
UnitCell, 8, 92
Label, 88
Spin, 89 central region, 43, 47, 49, 61
MPSH::, 17 checking the geometry, 14
Calculate, 89 command line options, see options, 79
128 Index
command prompt, 3 force, 6, 13, 25, 28, 108

comments, 36, 84 tolerance, 26
complex contour, 46 unit, 85
conductance, 5, 16, 46, 57, 85, 112 FracRecip, 11, 85
differential, 112 fractional
unit, 85 length unit, 91
constrain atoms, 26 reciprocal length unit, 11, 85
continuation symbol, 114, 116 spin unit, 85
contour integration, 46
convergence, 24, 28, 38 gate electrode, 20
problems, 25 GGA-PBE, 27, 38, 96
current–voltage characteristics, 17
HBar (spin unit), 85
default parameters, 9, 86–104 heterogeneous electrodes, 51
defaults file, 105 HOMO, 25, 111
density functional theory (DFT), 2, 23, 38 I–V curve, 17
density of states (DOS), 17 initial guess, 18, 50
DZ basis set, 28, 94 input and output files, 80
DZDP basis set, 28, 95 interfaces, 59
DZP basis set, 23, 28, 95 InverseLengthUnit, 85
effective potential, 19 k-point sampling, 36, 39, 44, 54, 62
electrode, 11–12, 42, 47, 49, 60 keywords, 86–104
alignment, 44
heterogeneous, 51 lattice constant, 8, 11
NetCDF file, 12, 14, 43, 61 lattice vectors, 8, 36
Electrode::, 12, 42 LDA-PZ, 23, 27, 31, 38, 96
AtomList, 12, 42, 92 LengthUnit, 84
ElectronTemperature, 61, 92 logical variables, 84
FermiTemperature, 72 LSDA-PZ, 27, 61, 96
NumKPoints, 12, 42, 48, 92, 93 LUMO, 25, 111
NumRep, 47, 52, 93
UnitCell, 42, 93 mesh cutoff, 28
electron temperature Method::Type, 93
bulk, 39 miswrapped electrode atoms, 45
molecule, 25 mixing, 24, 25, 98
energy molecular orbital, 7
Fermi, see Fermi level molecular systems, 6–7, 23–30
total, see total energy molecular unit cell, 29
unit, 85 Molecule::
zero-level, see zero-level, energy AtomList, 6, 94
environment variables, 81 ElectronTemperature, 25, 94
equivalent atoms, 44, 49 ElectrostaticPaddingFactor, 29, 72, 94
exchange–correlation functional, 2, 23, 27 FermiTemperature, 72
TotalSpin, 31, 94
Fermi level, 10, 39, 47, 103, 108, 111 Monkhorst–Pack, 36, 91, 92
Index 129
MPSH, 16, 21, 111 parallel calculations, 55, 59

Mulliken population, 63, 89, 106, 109 path, 81
multigrid method, 51 periodic systems, 7–11, 35–40
Poisson equation, 51
NEGF, 46 pseudopotential, 77, 108
NetCDF file, 7, 10, 14, 17 Pulay mixing, 24
NumOrb::, 93
BasisSet:: reciprocal lattice vectors, 8, 11, 36
Charge, 76, 95 Relaxation::, 25
DeltaRinn, 76, 95 ForceTolerance, 13, 26, 97
EnergyShift, 76, 95 MaxDisplacement, 26, 97
PseudoPotentialDirectory, 72, 77, 95 MaxSteps, 13, 26, 97
PseudoPotentialFile::XX, 77, 96 Type, 97
Size, 9, 13, 23, 28, 38, 72, 94 relaxing the geometry, 6, 13, 19, 25–26
SplitNorm, 76, 95 Results::, 105–112
V0, 76, 95 BandLines, 11
BasisType, 72 TransmissionSpectrum, 18
MeshCutoff, 9, 28, 96
RadialSampling::DR, 96 scaled coordinates, 36, 91
Tables:: scattering region, 41
Cutoff, 96 SCF::, 24
NumPoints, 97 Algorithm, 25, 97
XCFunctional, 23, 31, 61, 96 Beta, 72
Coordinate, 24, 25, 72, 97
options Criterion, 98
-a, 7, 15, 19, 79 DiagonalMixingParameter, 24, 25, 55,
-c, 14, 44, 79 72, 98
-h, 79 HistorySteps, 25, 55, 62, 72, 98
-i, 18, 50, 80 InitialBulkRun, 55
-o, 14, 79 InitialSpinPolarization, 31, 61, 98
-q, 80 MaxSteps, 25, 99
-v, 80 MixH, 72
--analyze, 7, 15, 19, 79 NumStepsReset, 72
--check, 14, 44, 79 Tolerance, 25, 97, 106
--force, 6, 14, 79 TwoProbe::InitialBulkRun, 99
--help, 79 screening, 41
--init-from, 18, 50, 80 scripts, 10, 113–119
--netcdf-file, 14, 79 self-consistency, 24
--param, 10, 28, 79 tolerance, 25, 98
--quiet, 80 SGGA-PBE, 27, 31, 96
--verbose, 80 Simulation::
--version, 79 ConstrainAtoms, 26, 100
--xyz, 14, 44, 80 Type::, 100
output file, 6, 9, 14, 105–112 Relaxation, 6, 13, 20, 26, 100
overwriting output files, 6 SingleConfiguration, 100
130 Index
spatially resolved density of states, 17 EquivalentCentralAtom, 44, 101

spin polarization, 30–33, 59 EquivalentElectrodeAtom, 44, 102
spin unit, 85 NetCDFfile, 14, 44, 73, 102
stacking, 47 Voltage, 50, 102
supercell, 48 RightSurface::NumberOfAtoms, 14, 16,
supercell, electrodes, 12 20, 102, 112
surface atoms, 20 UseMultigridForElectrostatics, 44, 51,
System::Type, 100 102
Bulk, 8, 36, 38
Electrode, 12, 42, 61 unit cell, 8
Molecule, 6 bulk, 8, 36
TwoProbe, 14, 44 electrode, 12, 42
SZ basis set, 9, 28, 94 molecule, 29
two-probe, 14, 44
SZP basis set, 28, 94
UnitCell::
TimeUnit, 85 LatticeConstant, 8, 10, 103
tolerance X0, 103
force, 26 Y0, 104
self-consistency, 25, 98 Z0, 104
total energy, 6, 9–11, 28, 29, 38, 40, 63,
vacuum level, 30, 111
100, 107
voltage drop, 19
transmission eigenchannels, 90, 112
transmission spectrum, 15, 45, 62 xyz file, 14, 44
two-probe systems, 11–21, 41–57, 59–67
TwoProbe:: zero-level, energy
CentralAtomList, 14, 16, 44, 100 bulk, 111
Contour:: molecule, 30, 111
Circle::NumPoints, 103 two-probe, 111
DeltaEMin, 46, 47, 55, 103
FermiLine::NumPoints, 103
KTPoles::NumPoints, 103
NumPoints, 46
RealAxis::Eta, 47, 103
RealAxis::PointDensity, 47, 103
LeftElectrode::
ATKfile, 72, 101
EquivalentCentralAtom, 44, 101
EquivalentElectrodeAtom, 42, 44, 101
NetCDFfile, 14, 44, 72, 101
Voltage, 17, 20, 50, 101
LeftSurface::NumberOfAtoms, 14, 16,
20, 101, 112
Molecule::GateVoltage, 20, 101
RightElectrode::
ATKfile, 73, 101
Atomistix A/S | Juliane Maries Vej 30 | DK–2100 Copenhagen, Denmark
Tel: +45 3532 0630 | Fax: +45 3532 0635

Tutorial and Reference Guide: Tomistix OOL IT

Uploaded by

Copyright:

Available Formats

Tutorial and Reference Guide: Tomistix OOL IT

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Tutorial and Reference Guide: Tomistix OOL IT

Uploaded by

Copyright:

Available Formats

ATOMISTIX TOOLKIT

Tutorial and Reference Guide

Atomistixr and Virtual NanoLabr are registered trademarks of Atomistix A/S.

Copyright c 2003–2006 Atomistix A/S.

Revision date: April 24, 2006 (2.0.3)

 ATK version 2.0, Atomistix A/S (www.atomistix.com).

 J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and

 J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 245407 (2001).

Remember to visit our website http://www.atomistix.com for updates to this manual

Atomistix, April 2006.

1.1 About This Tutorial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1 Molecular Orbitals of an H2 Molecule . . . . . . . . . . . . . . . . . . . 6

2.2 Band Structure of a Lithium Chain . . . . . . . . . . . . . . . . . . . . . 7

2.3 Setting up the Two-Probe System . . . . . . . . . . . . . . . . . . . . . . 11

2.3.1 Electrode Calculation . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3.2 Determining the Li–H Bond Length . . . . . . . . . . . . . . . . 12

2.3.3 Placing H2 Between the Lithium Chains . . . . . . . . . . . . . . 13

2.3.4 Self-Consistent Calculation at Zero Bias . . . . . . . . . . . . . . 14

2.4 Analyzing the Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5.2 Density of States . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.6 Current–Voltage Characteristics . . . . . . . . . . . . . . . . . . . . . . 17

2.6.1 Voltage Drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.7 Relaxation of a Two-Probe System . . . . . . . . . . . . . . . . . . . . . 19

2.8 Including a Gate Electrode . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.1 Water Molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.2 The Self-Consistent Loop . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.2.1 How To Cure Convergence Problems . . . . . . . . . . . . . . . 25

3.4 Changing the Accuracy of the Calculation . . . . . . . . . . . . . . . . . 27

3.5 Molecular Unit Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.6 Spin-Polarized Calculations . . . . . . . . . . . . . . . . . . . . . . . . 30

4.1 Band Structure of a Si Crystal . . . . . . . . . . . . . . . . . . . . . . . . 35

4.2 Convergence of an Al Wire . . . . . . . . . . . . . . . . . . . . . . . . . 38

5.1 Geometry of an Ideal Lithium BCC Lattice . . . . . . . . . . . . . . . . 41

5.1.1 Lithium Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.1.2 Two-Probe Geometry . . . . . . . . . . . . . . . . . . . . . . . . 43

5.1.3 Miswrapped Electrode Atoms . . . . . . . . . . . . . . . . . . . 45

5.1.4 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . 45

5.1.5 Energy Integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.2 Benzene Ring Coupled to Lithium Electrodes . . . . . . . . . . . . . . . 46

5.2.2 Central Region . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

5.2.3 Finite Bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.3 Heterogeneous Electrodes . . . . . . . . . . . . . . . . . . . . . . . . . 51

6 Spin-Polarization and k-point Sampling of Two-Probe Systems 59

6.1 Transmission Through a Fe–MgO–Fe Interface . . . . . . . . . . . . . . 59

6.1.2 Fe-MgO-Fe Geometry . . . . . . . . . . . . . . . . . . . . . . . 61

6.1.3 Transmission Spectrum . . . . . . . . . . . . . . . . . . . . . . . 62

6.2 k-point Sampling in Two-Probe Systems . . . . . . . . . . . . . . . . . . 62

6.3 k-point Sampling of the Transmission Spectrum . . . . . . . . . . . . . . 65

6.3.1 Plotting the k-dependent Transmission . . . . . . . . . . . . . . 65

Appendix A ATK vs. TranSIESTA-C 69

A.1 Changes to Keywords . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

A.1.1 Keywords No Longer Available . . . . . . . . . . . . . . . . . . 71

A.1.2 Keywords With New Names . . . . . . . . . . . . . . . . . . . . 72

A.1.3 Keywords With New Names and Syntax . . . . . . . . . . . . . 72

Appendix B Basis Sets 75

B.1 Basis Set Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

B.1.1 Single Zeta Orbitals . . . . . . . . . . . . . . . . . . . . . . . . . 75

B.1.2 Double Zeta Orbitals . . . . . . . . . . . . . . . . . . . . . . . . 76

ATK version 2.0, Atomistix A/S (www.atomistix.com).

J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and

J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 245407 (2001).