Home

Probe Particle Scanning Tunnelling Microscopy (PPSTM) simulation package:

Overview:

A Python3/C++ software package for simulating a couple of Scanning Tunnelling Microscopy - STM techniques based on results of electronic-structure computations (Density Functional Theory - DFT) using local basis set (LCAO). New users should have some experience with DFT computations and some overview of STM (e.g. 2nd chapter here)

Compilation and Installation:

Look at a dedicated wikipage for details of installation: installation wiki

GUI for a more user-friendly experience:

And follow the instructions from the full description of GUI

Developers:

Mainly developed at: Nanosurf Group, Institute of Physics of the Czech Academy of Sciences, Czech Republic and SIN group, Aalto University, Finland; with help from: Nanotech@surface Laboratory, EMPA, Switzerland. Contacts

Literature:

The structure of the code and how it works is explained mainly in the literature:

Main principle of STM and dI/dV calculations: Ondrej Krejčí, Prokop Hapala, Martin Ondráček, and Pavel Jelínek, Principles and simulations of high-resolution STM imaging with a flexible tip apex, Phys. Rev. B 95, 045407 – Published 6 January 2017 https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.045407
Main principle of IETS calculations and usage of d-sample orbitals: Bruno de la Torre, Martin Švec, Giuseppe Foti, Ondřej Krejčí, Prokop Hapala, Aran Garcia-Lekue, Thomas Frederiksen, Radek Zbořil, Andres Arnau, Héctor Vázquez, and Pavel Jelínek, Submolecular Resolution by Variation of the Inelastic Electron Tunneling Spectroscopy Amplitude and its Relation to the AFM/STM Signal, Phys. Rev. Lett. 119, 166001 – Published 16 October 2017 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.166001
More detailed description in a thesis of Ondřej Krejčí: DFT simulations of interaction of organic molecules with oriented surfaces, chapters 4 and 8. These chapters also contain the original papers above.
For simulations of images obtained with flexible tip apexes, please read about the Probe Particle Model (PP-AFM), whose calculations are necessary for these types of simulations: https://github.com/ProkopHapala/ProbeParticleModel/wiki and https://nanosurf.fzu.cz/wiki/doku.php?id=probe_particle_model (To be updated later)

Preparing Inputs with DFT codes:

Suggestion: It is better to run a geometry optimisation in a first run and then create the inputs for STM calculations, and possibly AFM pre-calculation, in a separate single-point calculation using the optimised geometry.

How to create the input files for the STM part of the Probe Particle is written here. You can also look at the ASE useful scripts for writing/creation the input geometry files.

This is also true for calculations with the 'fixed tip'.

For scanning with 'flexible tip apices' and storing in xsf and npy format, you must have PPAFM(>=0.2.0a3) installed. You can easily do it by running: pip install ppafm>=0.2.0a3.

For additional pre-calculations for the 'flexible tip' simulation and on the best strategy for the simulations, look here.

PPSTM_simple.py

Overview:

Fully automatic (closed) Python script for STM and/or dIdV calculations with stiff single atom or flexible tip-apex for calculations of systems on metallic/semiconductor substrates for sample biases between (-2,+2) volts:

You can copy this script to a directory with your (DFT and/or PP-AFM) pre-calculations, properly set ppstm_path to your PPSTM code (DO NOT USE $HOME for the path; both absolute and relative paths work), and you can run the script in the folder. Later you will have a possibility to see what options you used for the calculations.

The first part of the script contains variables that you can change as you prefer. There is a brief explanation for each variable at the end of the line in comments. The script should be self-explanatory.

All the Voltages are defined with respect to the sample (positive ~ empty states, negative ~ filled states) and the Fermi Level.

DFT codes inputs:

The script can handle DFT inputs pre-calculated from Fireball, CP2K, FHI-AIMS, and GPAW (in LCAO mode) codes. How to prepare those inputs is explained here. The script was tested for the first three DFT codes (feel free to test GPAW, but it has known problems for handling d orbitals on the sample and if any Hydrogens are not at the end of contributing atoms).

For CP2K and FHI-AIMS, spin-polarized calculations can be handled as well. You can choose 'alpha' = 'up' spin, 'beta' = 'down' spin, or 'both' spins contributing to your current.

Types of calculations/scans:

didv - single dI/dV (conductance) calculation at the energy (relative to the Fermi Level) given by V.

USE 'V-scan' instead of 'STM': 'STM' is calculated via rectangular integration = SUM of single dIdV for a single voltages between required voltage and zero. At the end, it gives you data for a single voltage V-scan: it calculates all dIdVs between Vmin and Vmax ( Vmin<=0<=Vmax ), then calculates STMs, and finally gives you all of them. You'll therefore have everything between Vmin and Vmax in a single run. If an experimentalist asks you about an image at 1.5 V, I would perform a V-scan between 0 and 2 V and then look at the image and compare. It is possible that either gap, the Fermi level, or some state is at the wrong energy from your DFT calculations, and therefore a better STM image can be found at a slightly different voltage. (Finally, if you want to go over +-2V, you have to adjust your cut_min and cut_max accordingly, and especially in positive voltage pray that your God of DFT makes your electronic structure of empty states far above the Fermi Level reasonable.)

For molecules: New feature 'states': it will make dI/dV calculations just for such Voltages that correspond to the eigenenergies (states = ei) of the sample, which are within V <= ei <= Vmax. You do not need to specify those eigen energies by yourself; the script finds them by itself. Beware that all the states are still taken into account for the dI/dV calculations. If you want to use just 1 state, choose very low eta (e.g. 0.0000001).

Periodic Boundary Conditions:

Various Periodic Boundary Conditions (PBC) can be used, but beware:

i) the script cannot use multiple k-points (it can read only Gamma k-point data). Therefore, the used supercell has to be big enough so the Gamma point will describe the electronic state of your sample properly.

ii) The script does not fold the scanning points into the original cell. Therefore, you can (sometimes) find the 'ends' of your geometry with 'boundary effects' taking place there. If the electronic structure is right (!), those effects cannot penetrate the inner part of the original super cell, but they affect a couple of Å near the end of the geometry.

iii) If the electronic structure is OK, you can: a) enlarge the number of multiplied cells and use the inner part only (possible with the PPSTM_simply.py script). Or b) come up with a procedure that will fold the positions into the original supercell (not written yet, one of the TODO goals).

Meaning of sample orbitals:

Orbitals of the sample are not directly connected with the used LCAO basis-set of your DFT calculations. Even though you can use very big basis set for any element, the PPSTM code uses only valence orbitals. Therefore you can use only sp sample orbitals, if your top-most atoms (those which are contributing to the current) are up to 3rd period of elements (even though i.e. you are using spd basis-set orbitals in Fireball.) Since d orbitals of the sample decay to the vacuum faster than sp orbitals, they are 'dampened' using constant rescaling of 0.2. (The factor is taken from comparison of FHI-AIMS s and d orbitals far away. The rescaling is not ideal, but different DFT codes use different decays for various elements. Therefore, we found no other simple general rule within the philosophy of this code.) Usage of sp orbitals is 2x faster than spd, and more tip-atom orbitals are programmed. You can also investigate the difference between sp and spd via comparing tip_orb = s' dIdV scans for a single height in your calculations.

Advanced features:

cut_atoms option can be used for (considerably) speeding up the calculations (especially for thick slabs and high pbc options). In most cases, you can cut out all the atoms that are 3 Å below the top-most (current-important) atom. (By current important, we mean anything that has considerable states on it. I.e., hydrogens on the flat organic structures have probably much lower contributions than C atoms, then the C atoms count). Or sometimes hydrocarbons on silicon can have much lower contributions than Si adatoms, then the Si adatoms count.)

If you have high-electro negativity differences in between atoms contributing to the current always check Hartree potential above these atoms. If big changes of the potential above your top-most atoms appear, try to estimate lower_atoms and_lower_coefs_ (rescaling for highly electronegative atoms). E.g., for O in the carbonyl group of the TOAT molecule, the oxygens' lower_coefs were set to 0.5 ( PRB 95, 045407 (2017) ).

Outputs:

Outputs can be in:

grayscale png images -- 2D for each height.
WSxM xyz files -- 2D for each height.
3D XSF files with full stack of data.
3D machine readable npy files with full stack of data.

Outputs can be for different heights, voltages, scanning techniques (STM or dIdV), and tip-atom orbitals:

Tip model and its orbitals:

Tip-atom orbitals are: s, px and py, pz, dz2, dxz and dyz. The last two options work only for sp orbitals of the sample, while the first three work with both sp and spd orbitals of the sample.

s tip-atom orbital works almost as Tersoff-Hamann theory (PRL 50, 1998 (1983) ) -- it goes beyond this theory, since it takes into account interference effects (as e.g. in ACS Nano 9 (9), pp 9180–9187 (2015) ).

We do not calculate the current for any DFT-calculated tip model. We rather opt for calculations with different tip-atom orbitals and comparison with experiment. For calculations with CO-tip which can be anything between s and pxy tip-atom orbital, depending on sample-bias, tip-base element and structure (from PRL 107, 086101 (2011) -- 50/50 s/pxy -- to PRB 95, 045407 (2017) -- 100% pxy tip) you can try some options like 'spxy' (50/50), 'CO' (13/87 -- PRL 119, 166001 (2017) ), '10spxy' (10/90), '5spxy' (5/95), but better is to calculate s and pxy in two separated runs and save it as npy or xsf. Then you can use the SUM_Orb_Contrib.py script to print you any configurations between s and pxy tip-atom orbitals. This is possible, since contributions from different tip-atom orbitals can be simply summed up.

All the 'fixed' scans are constant-height scans. Those should be done ~ 4 - 6 Å above the highest atom of the sample. To run a single-height scan, use: [zmin, zmax>zmax+dz, dz]. For the constant-current scans, there is no inner procedure. You can try to run scan for at least 2 Å heights with dz between 0.1 and 0.2, save it as xsf and plot it in VESTA or some similar software. But before you do this try to run a single constant-height scan if your simulations are similar to the experiment. Unless there is some crazy behavior with contrast reverse depending on height, there is no probable way that constant-current simulation would differ a lot from the constant-height one. Also, since I or dIdV ~ Exp(-beta z) you can try to plot the images in a log-scale (Not programmed in PPSTM_simple.py, you can use WSXM for that) to emulate the constant-current scan.

Atoms model plotting:

You can plot the position of atoms into the png figures and into the xsf files, if you have geometry, which you want to plot in an input_plot.xyz file. It can be 'bas' geometry, or it can be 'xyz' geometry with a single empty line bellow the amount of atoms in the file.

GUI for PPSTM code

USAGE

After cloning the repository and compilation:

Go to the PPSTM directory
Run > python GUI.py

The following interface should appear:

EXAMPLE

In the PPSTM/examples/ directory, you can find folders with example data, e.g., 4N-coronene. Inside the examples folder, you will find the already saved GUI options for this example: 4nCoroneneOpt.txt

Now, having the GUI in front of you, click Load options, go to the examples' directory and choose a .txt file with options. Everything is pre-set and should run without bugs
After loading the options, click the Import button. Wait for the terminal to say energies prepared, coefficients read
Now click Run

You've successfully run the simulation and plotted the image.

Below you will find a more detailed description of each GUI section with the instructions.

INPUT PART

Give required input parameters. All the fields in this section are necessary.

codes ~ You can choose from Fireball, FHI-AIMS, CP2K, and GPAW
sample orbital ~ either sp or spd
spin ~ none for spin-polarized calculations; up, down, or both for spin-polarized calculations
pbc ~ Affects boundary condition. The code is non-periodic from the nature of its design: 00 - calculates current from the original cell only - can be used for a single molecule, or if we want small part of big system; 11 calculated the nature from 3x3 super-cell centred around the original cell; 22 calculated the nature from 3x3 super-cell centred around the original cell - needed if the original cell is very small and we want to calculate big image.
lvs ~ 2D lattice vectors, works as 2D array or a file e.g. input.lvs
Give Input files path which is the path to your directory with necessary data: DFT pre-calculated eigen-energies and eigen-vectors and geometry
Give the name of geometry file, together with extension, e.g. input.xyz.
Give CP2K/GPAW ~ name used in CP2K or GPAW calculations
cut_at ~ number of atoms that contribute to the tunneling (-1 = all atoms)
lower_atoms ~ works as 1D array of integers; [] = None means that no atoms has lowered hopping; Be aware that python numbering occurs here: [0] -- means lowering of the 1st atom; [0,1,2,3] -- lowering of 1st 4 atoms
lower_coefs ~ works as 1D array of floats; [] = None means that no lowering of the hoppings; [0.5] -- lowering of the 1st atom hopping to 0.5; [0.5,0.5,0.5,0.5] -- lowering of 1st 4 atoms to 0.5

Now click Import ~ without importing, you won't be able to run the program.

If everything went correctly, the last line of the output should print: energies prepared, coefficients read

RUNNING OPTIONS PART

1st line: format, K, Q

Those parameters are used for a relaxed scan (i.e., tip type = relaxed) for which you need to have pre-calculated tip positions in Kx.xxQx.xx/PPpos_*.xsf or .npy files.

Format ~ choose if you need tip-positions as .xsf or .npy files

K ~ radial stiffness

Q ~ charge

Both Q and K are used for PPAFM calculations.
The following 3 lines of parameters are necessary for creating a rectangular cuboid grid. 0, 0, 0 is the beginning of the Euclidean system. The height of the scan should be about 4-6 Å above the highest atom of your system
Next 3 lines consist of parameters for the simulation

Scan type ~ there are 3 options: didv, states, and v-scan. The default v-scan is preferred, since it allows you to see all didv and STM maps for the <V,0> interval. states are useful for molecular systems - give all dIdV maps at energies of all eigenstates between Vmin and Vmax . didv calculates a single dIdV map at the energy of V.

Eta ~ Width of the Lorentzian function for energy smearing. Depending on the system, it can reach various values: for single molecular orbital very low number – 1e-6 eV; For standard slabs – 0.05-0.1 eV; For low layered or small slabs – 0.3-0.5 eV

WF_decay ~ Describes how the work function (barrier and therefore tunneling decay) is changing with increasing/decreasing voltage for each dI/dV step; 1.0 – the change of the tunneling decay scales with the voltage (at V=1.0 is WF=4.0 eV); 0.0 – no change with the voltage

V(min), Vmax ~ aplied sample bias = (energy vs. the Fermi Level in eV)

dv ~ It's a step between single dIdV calculations in case of V-scan. dV~eta

tip orb s ~ contribution to the tunneling to/from s orbital on the PP

tip orb pxy ~ contribution to the tunneling to/from non-tilting px and py orbitals on the PP (px contribution = py contribution = tip orb pxy /2)

OMP_NUM_THREADS ~ number of cpu cores for OMP paralelization: ncpu = 1 -- serial compilation & calculations; if ncpu > 1, then OMP paralel recompilation is used and C++ calculations are running on more cores

Tip type ~ fixed/relaxed, corresponds to fixed or relaxed scan
Now click Run

The main part of the PPSTM code will run, and the program will plot the first image. (By default it is: didv, height:0, voltage: 0. You can change that before clicking Run if you want)

PLOTTING PART

You can now change the plotted map with the following steps:

choose map type: STM / didv

STM is only available if the scan type was chosen to be v-scan
choose height and voltage index

Voltage index corresponds to different voltages if the scan type was chosen to be v-scan or states
Now click Show

You will see that the title of the plot changed to display your chosen map parameters

You can now save the image or the calculated data with the save image and save data buttons, respectively:

After clicking the button, in both cases the saving file window will pop up
You can save the image as a .png file
To save the data, you can choose a .xsf, .npy, or .xyz extension for it
By default, in both cases, the name of the saved file consists of data such as: map type, tip type, WF value, eta value, height, and voltage index

To save your chosen parameters (from both the input part and the running options part), you can click the save options button

The options will be saved as a .txt file with each parameter listed on a new line. You can write your comments into this file by making a new line and beginning a comment with #

Load options button lets you load those saved options back into the GUI. After clicking it, the file browser will pop up. You just need to select the file with your saved options and press Enter.

RELAXED SCAN

To run a relaxed scan, you need to have installed the PPAFM code with pip install ppafm>=0.2.0a3.

ADVANCED OPTIONS

If you need to change other parameters, you can inspect them inside pyPPSTM/guiMethods.py

Functions inside guiMethods.py:

imported() ~ code inside this function runs when you click the import button
newPPSTM_simple() ~ this corresponds to the Run button, and it calculates the main results

Scripting and Development:

Except for the 'closed' PPSTM_simple.py script and GUI, there is always a possibility for your own scripting. (But beware of necessary debugging):

You can look into a Si 7x7 example script to see how to calculated STM and dI/dV images with 'fixed tip'
An example how to run 'flexible tip' simulations with whole PP-AFM pre-calculations and the PP-STM calculations is in examples/4N-coronene/run_test.sh_. The script how to calculate the PP-dI/dV is in dIdV_test_4N-coronene.py and dIdV_test_4N-coronene.py_cp2k. For this kind of calculations a linking with PP-AFM code needs to be done: To read the positions of the relaxing PP the 'ppafm.io' are important. PP-AFM linking is also needed once you want to plot df together with STM.
To understand the input parameters of the inner dIdV and STM calculating procedure, look into input parameters of the running procedures.
To understand how the input files (from DFT codes) are read, look into reading parameters.

Usage examples:

Example of STM with rigid tip Si (111) 7×7 reconstruction: Si_7x7
Example of dI/dV scans above a spin-polarized CuPc molecule with a rigid tip: CuPc
Example of PPSTM simulations with flexible CO tip: 4N-coronene
Example of PP-IETS simulation with flexible CO tip: FePc_Au

TODO :

(but I will probably not have time to do it any time soon. Please send me a message if something would be interesting for you, I can try to speed it up.)