Table of Contents

Introduction

Probe particle model code is supposed to simulate AFM ( and to some extent STM and inelastic-STM) images obtained with tip modified by an atom or small molecule (such as Xe, CO, CH4, H2). In all cases the actual particle (atom or molecule) is replaced by a spherical model particle and classical potential (Lenard-Jones) is used for description of Pauli repulsion and Van der Waals attraction.

New code is written in C/Python and can operate in framework of Lennard-Jones forces as well as electrostatic forces, if necessary.

Older Fortran Version

Older documentation

New Version

The documentation is for “Master” branch, that can be downloaded from http://github.com/ProkopHapala/ProbeParticleModel/ or by running:

git clone http://github.com/ProkopHapala/ProbeParticleModel/

in your terminal.

It should be working for developers “dev” branch, too.

It can be downloaded from http://github.com/ProkopHapala/ProbeParticleModel/tree/dev or simply from terminal by running this commands:

git clone http://github.com/ProkopHapala/ProbeParticleModel/
cd ProbeParticleModel
git checkout dev

Software requirements

python-2.7; python-2.7-numpy(-1.8); gcc(-4.8); python-2.7-matplotlib(-1.3) - for plotting of figures;

Metacentrum
module add python27-modules-intel
module add scipy-0.17.1-py2.7.10

And don't forget to add:

import matplotlib as mpl;  mpl.use('Agg');

on the 3rd line of plot_results.py, if it is not allready in the script and if you want to use it on a cluster. This makes it run without Xserver (e.g. on supercomputer) # see http://stackoverflow.com/questions/4931376/generating-matplotlib-graphs-without-a-running-x-server

Hartree potential creation

DFT inputs

Inputs

Geometry

The geometry of the sample can be read from *.bas or *.xyz files, which has to have xyz extension! Here we show part of the geometry file from example attached to the model:

71

C       0.165   0.165   0.165   -0.0132352941
C       2.2987  1.39372 0.16643 -0.0132352941
C       4.4293  2.62743 0.16468 -0.0132352941
C       6.56208 3.85683 0.16431 -0.0132352941
C       8.69719 5.08888 0.16545 -0.0132352941
...

Which is in format:

Number of atoms optional line Symbol/or/Z-of-element x y z charge-optional

The Lennard-Jones (L-J) potential needed for calculations can be created by running:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/generateLJFF.py -i YOUR_INPUT_FILE.xyz

If charges are also in the input file add “-q” flag to create electrostatic field.

Optionally, the geometry can be also read from *.xsf or *.cube file, too. The command for creation of L-J force field then looks like:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/generateLJFF.py -i YOUR_INPUT_FILE.xsf

or

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/generateLJFF.py -i YOUR_INPUT_FILE.cube

x, y and z components of L-J forces are stored in LJFF_x.xsf, LJFF_y.xsf and LJFF_z.xsf files, respectively. These files that can be viewed e.g. via XCrySDen (http://www.xcrysden.org/) or VESTA (http://jp-minerals.org/vesta/en/).

Hartree potential

If an electrostatic Hartree potential is obtained from some DFT calculations, it can be read *.xsf or *.cube files. The electrostatic force field is created by running:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/generateElFF.py -i YOUR_INPUT_FILE

If default parameters are used, than you have monopole represented by an Gaussian cloud of charge with its FWHM of 0.7 Ǎ. The monopole can be changed to non-tilting dipoles or quadrupoles by adding flag: -t type, where type ∈ {s,px,py,pz,dx2,dy2,dz2,dxy,dxz,dyz}; s stands for monopole (default), p for dipoles, d for quadrupoles. The FWHM of the Gaussian cloud can be changed by adding flag: -s FWHM.

params.ini

This files contains all important information about the scan and informations for creation of important forcefields. Here we show an example of it:

probeType       8                               # atom type of ProbeParticle (to choose L-J potential ),e.g. 8 for CO, 54 for Xe  
tip            'dz2'                            # For calculations with electrostatics only - multipole of the PP {'dz2' is the most popular now fo CO}, charge cloud is not tilting  #
sigma           0.71                            # For calculations with electrostatics only - FWHM of the gaussian charge cloud {0.7 or 0.71 are standarts}  #
charge         -0.05                            # For calculations with electrostatics only: if 0.00 then ElFF is not even read - effective charge of probe particle [e] {for multipoles the real moment is q*sigma - dipole - or q*sigma**2 - quadrupole} {for CO 'dz2' we typically use -0.30 - -0.05}  #
stiffness       0.20 0.20 20.00                 # [N/m] harmonic spring potential (x,y,R) components, x,y is bending stiffness, R particle-tip bond-length stiffness, {for CO we typically use 0.24 0.24 20.00}
r0Probe         0.0 0.0  3.00                   # [Å] equilibirum position of probe particle (x,y,R) components, R is bond length {3.00 for CO mostly these days}, x,y introduce tip asymmetry
PBC             True                            # Periodic boundary conditions ? [ True/False ]
gridN           240 240 200                     # Grid division around each cell axis; Not necessary - if it is not here a 0.1 division is applied
gridA   23.10994667   0.000000      0.000000    # lattice vector of the cell; If geometry is read from .xsf or .cube
gridB   11.55497333  20.01380089    0.000000    # !!!! If geometry is read from *.xyz, but electrostatics from .xsf or .cube than gridN and gridA/B/C has
gridC    0.000000     0.000000    20.000000     # to be in agreement with the harteee potential. Also the shift vector in the cube file has to be 0.0 0.0 0.0 !!!! 
scanMin     00.0   00.0    10.0                 # start of scanning (x,y,z) (z should be something like: the highest atom of the sample + R + 2.5)
scanMax     30.0   22.0    15.0                 # end of scanning (x,y,z)
scanStep    0.1   0.1    0.10                   # steps of scan (dx, dy, dz)
Amplitude       1.0                             # [Å] oscillation amplitude for conversion Fz->df

If you want to make a scan for different probe, you have to change the probeType in params.ini and to recompute L-J forces.

The number of grid divisions in *.xsf files is enlarged by one in each direction. Therefore, gridN have to be numbers of cubicles in *.xsf file reduced by one, if geometry is read from *.xyz, but electrostatics from .xsf

Simulating AFM

With having L-J and electrostatic forces made by generating scripts, you can try to run an AFM scan via:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/relaxed_scan.py

which run the scan with charge Q and lateral stiffness K written in params.ini. Please note, that electrostatic forces are necessary only if Q ≠ 0.0. The result - Fz force acting on the tip – is saved in Q?.??K?.?? directory as an OutFz.xsf file. The df results for constant height scans can be plotted by:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/plot_results.py  --df

that plot results for oscillation written in params.ini. The results - df_???.png maps for each tip height - are in directory: Q?.??K?.??/Amp?.??

optionally WSxM files for each height can be written by putting flag - -WSxM behind the plotting command. Outputs are df_???.xyz files with WSxM head and x y df data. If - -save_df flag is applied the df data are stored in df.xsf file. Flag - -atoms used simultaneously with - -df puts positions of atoms of the sample saved in input_plot.xyz into df_???.png. The flag - -bonds ads into the maps also lines connecting close-by atoms.

If a flag - -pos is applied for both commands (scanning & plotting) than xy positions of the relaxing Probe Particle (PP) are shown in xy_???.png as a red dots, while the gray scale on the background maps represent z position of the PP (brighter - higher).

Tests

Examples of df simulation is already in the downloaded/cloned repository in folder examples/. You should try to run these examples before going to your own stuff, in order to see that the code is working on your machine.

Scans with different charge (Q), lateral stiffness (K) or oscillation amplitude (A)

A scan with different charge (Q) and/or lateral stiffness (K) than those written in params.ini can be calculated via running (be aware, that for Q ≠ 0.0, you need to have precalculated electrostatic forces):

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/relaxed_scan.py -q (Q) -k (K)

It will create a new folder Q?.??K?.??/ where will be stored the results from the new run. A df map for wanted oscillation amplitude (A) can be then created by:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/plot_results.py --df -q (Q) -k (K) -a (A)

The results will be saved in subfolder Amp?.??/.

Also scans over ranges of (Q),(K) and (A) can be done, via:

python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/relaxed_scan.py --krange min max nK  --qrange min max nQ
python PATH_TO_YOUR_PROBE_PARTICLE_MODEL/plot_results.py --df --krange min max nK  --qrange min max nQ --arange min max nA

References

Prokop Hapala, Georgy Kichin, Christian Wagner, F. Stefan Tautz, Ruslan Temirov, and Pavel Jelínek, Mechanism of high-resolution STM/AFM imaging with functionalized tips, Phys. Rev. B 90, 085421 – http://journals.aps.org/prb/abstract/10.1103/PhysRevB.90.085421

Prokop Hapala, Ruslan Temirov, F. Stefan Tautz, and Pavel Jelínek, Origin of High-Resolution IETS-STM Images of Organic Molecules with Functionalized Tips, Phys. Rev. Lett. 113, 226101 – http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.113.226101