Home

Probe Particle Model

This is implementation of efficient and simple model for simulation of High-resolution atomic force microscopy (AFM), scanning probe microscopy (STM) and inelastic tunnelling microscopy (IETS) images using classical forcefields. The main principle is illustrated in following image:

The schematics illustrates all the forces acting on the relaxing probe particle (PP - blue ball), which represents the last atom of the flexible tip apex. Those forces are: i) spring forces - representing the forces of the tip and keeping the PP bellow the last atom of the metallic tip-base (yellow ball); ii) Pauli repulsion and van der Waals attraction between the PP and fixed atoms of the substrate approximated by Lennard-Jones potential; iii) Electrostatic forces between sample (calculated via DFT) and (monopole/dipole/quadrupole) charged PP (if PP/flexible tip apex is charged). Detailed description of the forces and difference between the original paper and current implementation are described here

The code has a form of python module which imports C/C++ library using ctypes. This way it should be very friendly to any further development and integration to other projects.

Quick start

Installation & Compilation

UNIX/LINUX & partially MAC OS X users:

Simply clone the github repository using

> git clone https://github.com/ProkopHapala/ProbeParticleModel.git

The code depend only on Python 3 (tested with versions >=3.6) with standard numpy (>=1.16) and matplotlib (>=3.0) modules and g++ (>=5.4) compiler (Especially for MAC user: if you use different C/C++ compiler then you need to edit your copy of the makefile ). Otherwise it is completely self-contained.

MAC OS X additional comments: here.

Windows users: here.

The C++ part of code is automatically recompiled every time you import the library in python. This makes debugging easier since you are sure you work with up-to-date C++ code. On the other hand is may cause some problems for Windows or OS X. In case you need to switch it off please edit corresponding cpp_utils.make() procedure; more details about it can be found here.

Inputs

For the minimal input without any electrostatics only atomic positions are needed:

• They can be provided in a bas or xyz format in and *.xyzfile. !! Don't forget to clear the 2nd line of the xyz format !! . You can also add (e.g. Mulliken) charges as an 5th column to this file. See the forcefield section for more details. Also don't forget to specify the cell into the params.ini file.
• Another option is to use FHI-AIMS geometry format *.in. The cell can be read directly from the geometry file, if specified there.

If you have an possibility to calculate the electrostatic = Hartree potential from an DFT codes, this should give you more precise results. At the moment 4 codes are supported: Fireball, VASP, FHI-AIMS and CP2K. (Using other codes shouldn't be a big problem, but beware of different units. Please write us, if you are unsure, or you are running into problems.) The files can be used to provide everything - information about electrostatics, geometry and cell. Look at the forcefield section how to run the calculations using these.. With these 4 codes two input file types should be used:

• *.xsf -- from Fireball or VASP (use v2xsf to transfer LOCPOT into LOCPOT.xsf). Potential is in eV (same units as PP is using).
• *.cube -- from FHI-AIMS or CP2K. Potential in in Hartree. PP uses internal rescaling.
• Please do not mix these workflow-paths, otherwise you can encounter problem with units.

Examples

We have prepared a few examples which should run out-of-the-box under LINUX just by running run.sh bash script. This include

• very simple and fast calculation of Graphene with nitrogen defects using point-charge model of electrostatics
• A bit more sophisticated calculation of monolayer of PTCDA molecule on Ag111 with DFT Hartree potential as input. The Hartree potential in LOCPOT.xsf is downloaded automatically by wget from our hosting. (Usually doesn't work on OSX, it is better to download the file manually and update run.sh accordingly). There are two examples: first of them is using monopole or quadrupole charge on the probe particle. The differences between monopole and quadrupole charge is illustrated in figure from Supplementary Information of Nat. Commun. 7, 11560 (2016):
• IETS example recalculating the theoretical data without any electrostatics published in PRL 113, 226101 (2014).

Running calculations

As you can see from the content of the directories 2 input files are required:

1. params.ini which contains input parameters of simulation (dictionary of all parameters can be found in pyProbeParticle.common.params).
2. *.xyz which contains input atomistic geometry (first column atom type, followed by 3 columns of x,y,z cartesian coordinates of the atom, followed by optional column of atomic charges ) ... see example
   • In more sophisticated calculations you can use Hartree potential from DFT calculation instead of atomic charges, in this situation both the atomic geometry can be read from single *.xsf or *.cube file. Example of first of them can be seen while running the example in https://github.com/ProkopHapala/ProbeParticleModel/tree/master/examples/PTCDA_Hartree
   • Atomic geometry and parameters of the cell can be also read from geometry.in FHI-AIMS file.

As you can see in run.sh the calculation is composed of 3 independent steps:

1. Calculation of Lennard-Jones and Electrostatic forcefield over sample.

This is done by running script:

python generateLJFF.py -i input.xyz -q

if you use point-charge electrostatics model and geometry in input.xyz or by:

python generateElFF.py -i LOCPOT.xsf
python generateLJFF.py -i LOCPOT.xsf

if you have Hartree potential and geometry stored in LOCPOT.xsf file ( LOCPOT.xsf can be generated e.g. from VASP using v2xsf utility). Hartree potential from Fireball (fftpot.xsf) or generated from FHI-aims or CP2K codes (*.cube files). The computation of electrostatic force is mainly explained in Supplementary Material of PRL 113, 226101 (2014).

Output FFLJ_?.xsf is x,y,z component of Lenard-Jones force in units [eV/Angstroem] for posistion of Probe-Particle sample on grid defined by sampling of input LOCPOT.xsf or defined by gridA,gridB,gridC,gridN parameters in parameters.ini.

The FFel_?.xsf contains x,y,z component of Electrostatic filed in [V/Angstroem] (scaling by probe charge Q is done only in relaxed scan procedure)

2. Relaxed scan of tip over sample,

where position of probe particle is relaxed in order to minimize its total energy, i.e. to find balance of force between forces from tip and sample. The sample force is found by interpolation of total force-filed obtained Ftot = FFLJ + Q*FFel. Spring-potential from tip. This is done by executing script:

python ../../relaxed_scan.py -k 0.5 --qrange -0.05 0.0 2 --pos

where -k denotes tip stiffness in [N/m]; --qrange range of charges on probe particle; and --pos indicates that beside the force also relaxed positions of probe particle should be exported. The results for different stiffness and probe charge are stored to independent directories.

3. Plotting,

where we just read output OutFz.xsf and PPpos_?.xsf and generate set of images from it. Here we also convert of vertical force Fz to frequency shift df for given (peak-to-peak) oscillation amplitudes set by -a or --arange and given parameters of the cantilever kCantilever,f0Cantilever (default parameters are for standard Qplus sensor). This is done be executing script:

python ../../plot_results.py -k 0.5 --qrange -0.05 0.0 2 --arange 0.5 2.0 2 --pos --df 

Flag --df enables plotting of the df_???.png images with calculated df data, while --pos plots the positions of probe particle positions during the simulated scan into xy__???.png This is the only part of the procedure where dependence on matplotlib is required. Instead you can also use other visualisation tools:
- By adding --save_df flag you can save your dfresults in df.xsf file, which you can visualise e.g. via VESTA.
- Using --WSxM flag you will print each height slice into df_???.xyz file, readable e.g. by WSxM visualisation tools.

Further outputs

The z-component of the force acting on the tip can be plotted into Fz_???.png files by using --Fz flag. (Second electrostatics is not taken into account in the plotted force!)

The IETS signal (as calculated in PRL 113, 226101 (2014), that means using variation of CO frustrated translation peak energy approximation, see PRL 119, 166001 (2017) for more details) will be calculated and plotted via applying --iets=M,V,W flag, where: M - probe particle mass [a.u.]; V - bias offset [eV]; W - peak width [eV].

Both of these quantities can be printed into *.xyz files by means of --WSxM flag.

Optional flags:

• --cbar: adds a colorbar with absolute numbers to the *.png images.
• --atoms: adds a position of atoms saved in input_plot.xyz into the plotted *.png images.
• --bonds: adds also lines in between close-by atoms ...
• --no_int: allows for plotting of images without any interpolation between pixels

You can also speed up the calculations and spare approx. 1/2 of your disc memory by using --npy behind all the scripts. It will save all your intermediate data into "machine readable" *.npy files, which on the other hand are hard to be debugged.

There are also other flags, that are connected with other developments. For example see KPFM approximations here.

Note on parameters in params.in and double electrostatics:

We strongly recommend to keep all the parameters for the calculations (mainly multipole = tip, sigma, PBC and the second multipole = tip_base) in the params.ini instead of using flags.

Any line in the params.ini file can be commented by using # (ideally with space).

Multipole (dipole, quadrupole) charge is not multipole moment. Dipole moment can be calculated as charge * sigma [e*A]; quadrupole moment is charge * sigmaˆ2 [e*Aˆ2].

CO quadrupole charge vary from -0.05 to -0.30 (using sigma = 0.71 or 0.7) depending on experiment.

Stiffness of 0.24-0.25 N/m (stiffness 0.24 0.24 20.0) is mostly used for CO-tips (Nat. Commun. 7, 11560 (2016)).

PBC affects only calculations of L-J poetential. The code is always periodic from the nature of its design: PBC False means nPBC 0 0 0 ; PBC True means read nPBC parameter.

nPBC affects only calculations of L-J poetential -- how many cells around the orignal one are used for L-J force-field caclulation. The last number plays no-role. The total amount of cell is given by (nPBC(x).2+1).(nPBC(y).2+1) and the final geometry is always centred in the original geometry as given in the input file. nPBC 0 0 0 means only the original cell; nPBC 1 1 1 3x3 cells around the original one and so on. See the image bellow for graphical representation.

gridN can be assigned automatically (commented or missing in params.ini) with division of 0.1 A.

The scanning height is referred to the position of the last atom of the metal tip-base. The equilibrium position of the probe particle is then by r0Probe-z component lower. (4 A is default; e.g. in Nano Lett. 16, pp 1974–1980 (2016) and ACS Nano 12, 5274−5283 (2018) 3 A were used.)

If in IETS you get a loads of negative frequencies, try to change the scanning grid or gridN. There are always few negative energies, especially in really close distances.

Cell definitions

The code expects, that your system is in the cell starting in 0,0,0, given by the 3 lattice vectors - each of them is a row: 1st is gridA; 2nd gridB and the last one is gridC (normally 0. 0. height of cell ) - basically the same strategy as in ASE:

• For example if your system has centre of mass at 0,0,0 you can run into problems (sometimes problem with FHI-aims).
• Or if your geometry is bellow 0 in z then the atoms are not captured properly, regardless of stated PBC or nPBC.
• it is also good practise to leave enough space above the top-most atoms in the z direction: 15 Angstrom is pretty safe.
• which leads to another good practise to shift the system to a lower part of your cell.

If you are using hartree potential (xsf or cube files); Then you do not need to care about these - the points of the grid are exactly the same as in those file with potential.

• If you want better (more dense) spacing or larger cell, then you need to redo you DFT calculations with adjusted parameters

If not -- you have to specify: gridA , gridB and gridC by hand (default values do not make any sense):

• if you have system that is periodic in x and y (surface) than adjust gridA and gridB with the exact lattice vector; gridC needs to have enough space about the atoms (as mentioned above).

Double Electrostatics

Second electrostatics (ACS Nano 12, 5274−5283 (2018)) can be applied by keyword tip_base, which has two parameters: multipole and charge. (The original electrostatics is still handled by keywords tip and charge) I.e. the metal tip-base dipole, due to Smoluchowski effect, of 1.5 Debye (Nano Lett. 16, pp 1974–1980 (2016)) can be applied by using tip_base 'pz' 0.44. This line has to be in the params.ini during calculations of generateElFF.py, relaxed_scan.py and plot_results.py.

The actual metal tip-base charge is applied only in the last script (Fz -> df conversion only!) and does not have any effect on the relaxation or IETS!

Image illustrating electrostatic field around CO-tip - from left to right (1-3 and 5 are analytical models) 1) Electrostatic potential around quadrupole on oxygen. 2) Potential around dipole of the metallic tip-base. 3) Electrostatic potential of the whole system. 4) DFT calculated potential around CO-tip on metallic cluster, taken from Supplementary Information of Nat. Commun. 9, 122 (2018). 5) z-component of the electric field around the whole (CO-tip + metallic tip-base) system. Green point (at the corner of triangle) depicts the last atom of the metallic tip-base, while gray and orange points illustrates positions of carbon and oxygen, respectively. The orange circle represents van der Waals radius of the oxygen atom

Note: According to experience, the double electrostatics is important mainly for materials with strong polarization, while for (most of) molecules on metal substrates it has minor effect as you can see on PTCDA/Ag(111) simulations without/with (left/right) metal tip-base dipole and quadrupole: Q: -0.3/0.0/0.3 (top/middle/bottom), sigma 0.71.

Results Example

Example of Simulated (right) and experimental (left) High resolution AFM image of self assembled monolayer of PTCDA molecules on Ag(111) surface with Xe tip ( at liquid helium temperature with collibry senzor ). Positive monopole of 0.3 e was used for the simulations

References:

Main 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 – Published 19 August 2014

• 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 – Published 25 November 2014

OPLS force-field:

• W. L. Jorgensen and J. Tirado-Rives, The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin, J. Am. Chem. Soc. 110, 1657 (1988)

Using of multipoles (quadrupole):

• J. Peng, J. Guo, P. Hapala, D. Cao, R. Ma, B. Cheng, L. Xu, M. Ondraćěk, P. Jelínek, E. Wang, Y. Jiang, Weakly Perturbative Imaging of Interfacial Water with Submolecular Resolution by Atomic Force Microscopy. Nat. Commun. 2018, 9 (1), 122

Double electrostatics:

• Michael Ellner, Niko Pavliček, Pablo Pou†, Bruno Schuler, Nikolaj Moll, Gerhard Meyer, Leo Gross, and Rubén Peréz, The Electric Field of CO Tips and Its Relevance for Atomic Force Microscopy, Nano Lett., 2016, 16 (3), pp 1974–1980

• Fabian Schulz, Juha Ritala, Ondrej Krejčí, Ari Paavo Seitsonen, Adam S. Foster, and Peter Liljeroth, Elemental Identification by Combining Atomic Force Microscopy and Kelvin Probe Force Microscopy, ACS Nano 2018, 12, 5274−5283

How it works

The code consider basically tree types of interaction between tip ( PP particle) and sample.

• Pauli-repulsion ( currently approximated by repulsive part of Lenard-Jones potential r^-12 )
• van der Waals attraction ( currently approximated by attractive part of Lenard-Jones potential r^-6 )
• electrostatic ( currently can be computed as coulomb pairwise interaction between PP and point charges in centers of atoms of sample, or by reading electrostatic force-field obtained by derivative of sample Hartree potential as described in supplementary of PRL 113, 226101 (2014).

The computation of images is divided into two parts:

• precompute vector Forcefield ( Fx(x,y,z), Fy(x,y,z), Fz(x,y,z) ) over sample and store it on a 3D-grid. ( see getLenardJonesFF and getCoulombFF functions for more details ). After individial components of forcefiled are sampled ( i.e. Lenard-Jones and electrostatic ) they are summed up to form one single total effective forcefield in which Probe-Particle moves.
• relax Probe-Particle attached to the tip under influence of the total effective forcefield. The relaxation of Probe-Particle is done using "Fast Inertial Realxation Engine" (FIRE) algorithm ( see.
  PRL 97, 170201 (2006) implemented in FIRE:move() function. In each step of relaxation the forcefield on the grid is interpolated in interpolate3DvecWrap function, considering periodic boundary condition. From python the relaxation is called as relaxTipStroke function providing 1D-array of tip positions, and obtaining back 1D-array of Probe-Particle position after the relaxation and of force between tip ProbePartcile and sample at that position. The lateral scan in (x,y) to obtained stack of images is done by calling relaxTipStroke at different (x,y) position of tip, where each call of relaxTipStroke does one approach along z-direction. This correspond

Why it is splited like this ?

This splitting of computation ( first sampling of forcefiled on the grid, and than relaxation using interpolation) has several advantages over straightforward computation of interaction on the fly during relaxation process.

• It is faster - if there is ~100 atoms of sample, summing over each pairwise Lennard-Jones interactions, is much slower, than just interpolating the forcefield from the grid. Because force evaluation is done several-times for each voxel of the 3D scanning grid, it is more efficient to precompute it just once for each voxel, and then interpolate.
• It is more flexible and general - Decoupling of relaxation process and computation of forcefiled allows us to plug-in any forcefield. The original motivation was to used electrostatic forcefield obtained from Hartree potential from DFT calculation. However, it is not limited to that. We can plug in e.g. a derivative of potential energy of probe peraticle (i.e. Xe or CO ) obtained by scanning it in DFT code, in similar way as Guo did in JPC C 119, 3, 1483-1488 (2015). The only limitation here is computational cost of obtaining such potential from ab-initio calculation.