Nanosurf Lab

The code is written in Fortran 90 using object oriented features and custom types. In some parts OpenMP paralelization may be strightforward and efficient, but it is not implemented yet.

The source codes ( .zip)

1. Example CoPc -

older version without electrostatics does not need big .xsf files as input) … the code is not exactly described in this tutorial but is very similar (use pairwise versions of subroutiens instead of grid interpolation)

plug-and-play example with older code - just set proper path in SHTM.sh and run it

1. T_atoms.f90 … definition of 'atomType' ( i.e. chemical element with given set of parameters) and 'subSystem' (subsystem is piece of molecular geometry composed of atoms such as tip and surface geometry )
2. T_grid3D.f90 … definition of scalar rectangular uniform grid, implements subroutines for loading and writing .xsf files, interpolation, initialization …
3. T_grid3Dvec.f90 … the same for vector grids used to describe force fields ( such as electrostatic force field )

1. G_globals.f90 … there is defined a lot of global variables. Some of these could be possibly made local in future.

1. relax.f90 … implements:
   1. main relaxation MD loop 'relax()' and 'relaxGrid()', There are evaluated of forces ( either pairwise terms by analytic formula or by grid interpolation ). The 'relaxGrid()' is used in final grid-accelerated version. The relax() is there for historical reasons and for debuging.
   2. move_simple(), move_FIRE(), move_FIRE_true() subroutines which implements one step of relaxation algorithm ( are called from relaxation loops in each step ). Default and most effective is move_FIRE_true() which implements FIRE minimization algorithm ( see Bitzek, E., Koskinen, P., Gähler, F., Moseler, M. & Gumbsch, P. Structural relaxation made simple. Phys. Rev. Lett. 97, 170201 (2006). )
2. forceFields.f90 … implements:
   1. evaluation of functional form of various emperical forcefield evaluation like 'getFF_LJ()', 'getFF_Morse()', 'getFF_Harmonic()'
   2. evaluation of functional form of tight binding hopping elements (with s-s and p-p orbitals ) used for evaluation of STM images
3. sampling.f90 … 'sampleSurf()' do sampling of pairwise L-J potential from surface atoms for discrete grid of fixed Probe Particle position and store it into a grid3Dvec type variable for further use during actual simulation ( with relaxation).
4. dynmat.f90 … assembling and solving 3×3 dynamical matrix by sampling force field around given point. There are again 2 versions one for which evaluate surface potential as pairwise analytical terms directly (older version, debugging) and the newer one which use interpolation from precomputed grid (this is used by default).
5. eigenSolver.f90 … subroutine with stable analytical algorithm for solving eigenvalues of 3×3 matrix
6. readfiles.f90 … subroutines to read atom types from 'species.dat' and geometry from a .bas file ( tip.bas, surf.bas )

1. SHTM_3D.f90 … it is structured in following order:
   1. read input parameters from SHTM_3D.ini
   2. call readspecies() to load L-J parameters
   3. read geometry of tip and sample from files
   4. read 'poslist.ini' ( list of restarting position if Probe particle is stuck in unfavorable local minimum under pressure )
   5. creating grid force field. By defaul 'Felec_X.xsf' files are loaded, and L-J force field is sampled by 'sampleSurf()' (without relaxation - just for Probe Particle set to particular fixed grid points). This is often the most time-consuming step.
   6. initialization of output grids
   7. main loop for performing relaxed 3D scan. By default ( with onGrid = 1 ) relaxation is done using interpolation of precomputed grid force field. In this loop all properties ( such as tip deviation, Force, STM hopping and vibration eigenvalues ) are obtained. In case when onGrid = 1 is used, the main loop is often faster than evaluation of force field on grid.

• Warrning : Only rectangular supercells are curretly considered. If you use non-rectangular supercell results are unexpected.

1. Before computation using Probe Particle Code, it is necessary to prepare electrostatic force field for Hartree potential obtained from some DFT package ( like vasp, Fireball … ). For this purpose I prepared quite general script ( potToFF_final.py ) which reads file 'LOCPOT.xsf' and outputs files Fx.xsf,Fy.xsf,Fz.xsf containing x,y,z electorstatic force component acting on 3D-gaussian-shaped cloud of charge containing 1 electron in total. The scripts performs derivatives and convolution in Fourier space ( using Fast Fourier Transform ) for good performance and numerical robustness (instead of doing finite differences).
2. Prerequsities : python (tested with 2.7 ), numpy, pylab, matplotlib
3. Usage - put LOCPOT.xsf and potToFF_final.py into same directory and run it. It will automatically produce Fx.xsf,Fy.xsf,Fz.xsf files. There are no command line parameters, all settings should be editted inside the script.
4. The most relevant setting for user is sigma = 1.0 # [ Angstroem ] which determine width of gaussian cloud of charge. It is supposed to be set approximatively as a radius of an atom.
5. after generation the files Fx.xsf,Fy.xsf,Fz.xsf should be renamed to Felec_1.xsf, Felec_2.xsf, Felec_3.xsf and copied to work directory of the main simulation.

Python scripts (.zip) for conversion of Hartree DFT potential to grid force field. Detailed description of all scrips is beyond scope of this tutorial.

possible some scripts may need some python library from my this python path directory (wich I have on my computer). If so, just extract this file and copy the required *.py file into your working directory or python path. pokop_python_path.zip

There are described parameters of simulation

 0		                # 0..2 level of debugging output
 2           	                # output format 1=ppm 2=xsf 3=tab  ( input, output, Debug) ... mostly tested just with .xsf
 3 1		                # relaxation_method (0=fixed 1=damped_MD 2=FIRE_my  3=FIRE_true ) 3=FIRE_true is recommanded for best stability and performance   |   1 or 0  onGrid - if acceleration by precomputed grid should be used, always use 1 ( 0 is just for debugging )
 1  1000                      # Periodic Boundary Condition ( use 1 )     |  FgridMax ... limit maximum force evaluated on grid ( to protect divergence and numerical instability around centres of atoms  in eV/A )
 0.3 0.000001 0.1 0.0 10000   # dt [fs] | convergence criterium [eV/A] | damping | startKick (multiples of dt for first step of MD)| maxRelaxIter (if not converged until maxRelaxIter skip pixel )
 0.0000000          	        # kMorse [eV/A^2], if kMorse<=0 use Lenard-Jones instead ... currently not used
 1.0 1.0			# beta1 beta2 [1/A] - decay factor for STM simulation ( if just relative contrast of images does matter just use 1.0 )
 0.09 0.09 1.0    	        # kHarmonic(x,y,z), aditive harmonic potential of tip bending  in eV/A^2
 0.00 0.00 0.00               # lateral spring minimum ( to simulate tip asymmetry ) in [Angstroem]
 0.025			# displacelent for dynamical matrix [Angstroem]
 7 -0.4                 	# probe_atom_type ( refer to index in species.dat ) | probe_charge (in elementary charge )
 0.1 0.1   0.1		# scan step [Angstroem]    
 -2.0 -2.0   6.0	        # Scan RMin [Angstroem]
 23.0 23.0  12.0	        # Scan RMax [Angstroem]

There are described parameters of atom types. The first column is the index of atom type which should be referred to in surf.bas and tip.bas. Next is vdW radius in Angstroems (for two different atoms the equilibirumu distance is sum of this numbers). Third column is binding energy for two atoms of this type (for two different atoms the binding energy is geometrical average). 4-th column is just caption used when .xyz file is written (usually just for debugging ). behind '#' sign is just comment with reference of force field parameters used.

 8   # number of species (legend is at the end )
 1	1.487	0.0006808054	H	# HC abalone AMBER 2003  OPLS
 2	1.908	0.0037292524	C	# C abalone AMBER 2003  OPLS
 3	1.8240  0.0073719       N	# N abalone AMBER 2003  OPLS
 4	1.6612	0.009106314	O	# O abalone AMBER 2003  OPLS
 5	2.0000	0.030          	Co	# Cu test Xe-Cu binding 0.2 eV/A Silvestrelli, P. L., Ambrosetti, A., Grubisi, S. & Ancilotto, F. Adsorption of rare-gas atoms on Cu(111) and Pb(111) surfaces by van der Waals corrected density functional theory. Phys. Rev. B 85, 165405 (2012).
 6	1.6000	0.030          	Ag	# same as Cu
 7	1.6612	0.009106314	O	# O abalone AMBER 2003  OPLS
 8	2.0000	4.0000000000	Au	# AuTip
 # Z    R0[A]    E0[eV]        Name   #Comment 

contains description of positions of atoms on surface (the molecular layer perhaps). Note that first column should denote index of the atom type (as described in species.dat) not proton number of element.

  24 
  1     -0.39724     1.21664    -0.00000 
  1      0.39739    -1.21995     0.00000 
  2     -2.07709    -0.18301     0.00000 
  2     -3.54836    -0.13178     0.00000 
  1     -4.07478    -1.10172     0.00000 
  1     -4.19079     3.25336     0.00000 
  2     -3.60804     2.31915     0.00000 
  2     -2.23504     2.36291     0.00000 
  2     -4.24726     1.04736     0.00000 
  1     -5.35352     1.01238     0.00000 
  1      4.19973    -3.24228     0.00000 
  2      3.61338    -2.31031     0.00000 
  2      2.24052    -2.35982     0.00000 
  2      2.07269     0.18509     0.00000 
  2      3.54400     0.13997     0.00000 
  1      4.06600     1.11228     0.00000 
  1      5.35391    -0.99670     0.00000 
  2      4.24783    -1.03610     0.00000 
  3     -1.47862     1.17719     0.00000 
  3      1.47862    -1.17719     0.00000 
  4     -1.36960    -1.17719     0.00000 
  1     -1.66143     3.30156    -0.00000 
  4      1.36266     1.17729    -0.00000 
  1      1.67157    -3.30128     0.00000 

Description of geometry of tip apex. The code is made with a vision to enable usage of different tip geometry. However currently there is a limitation that probe particle position (and constraining lateral spring) always refer to the first atom in the list. Up to now I almost exclusively used just single atom tip apex. The first column again denote index of atomic type in 'species.dat'

 1
 8      0.00000      0.00000      0.000000

The file contains restarting positions where the code will try to initialize Probe Particle position if current postion has unfavorable energy ( for more details see do 'ipos = 1,npos' sub loop of main simulation loop ). This file is currently required but the functionality is controversial. May be removed in future.

First number denote number of positions, while the rest is list of tried positions described in spherical coordinates. First column is multiple of tip-to-probe equlibirum distance, second is polar angle and third is longitudinal angle.

 21
 1.0 0.0 0.0  # just under tip apex
 1.0 0.1 0.0
 1.0 0.1 0.1
 1.0 0.1 0.2
 1.0 0.1 0.3
 1.0 0.1 0.4
 1.0 0.1 0.5
 1.0 0.1 0.6
 1.0 0.1 0.7
 1.0 0.1 0.8
 1.0 0.1 0.9
 1.0 0.2 0.0
 1.0 0.2 0.1
 1.0 0.2 0.2
 1.0 0.2 0.3
 1.0 0.2 0.4
 1.0 0.2 0.5
 1.0 0.2 0.6
 1.0 0.2 0.7
 1.0 0.2 0.8
 1.0 0.2 0.9

1. FFgrid_1.xsf,FFgrid_2.xsf,FFgrid_2.xsf … the precomputed cartesian components of total forcefield ( Lenard-Jones + Electrostatics ) on grid (1=x,2=y,3=z) which is used for interpolation during the calculation
2. OutEig1.xsf, OutEig2.xsf, OutEig3.xsf … stiffnes eigenvalues in units of [ eV / A^2 ] form solution of dynamical matrix aster relaxation of probe-particle for each tip position
3. OutFx.xsf,OutFy.xsf,OutFz.xsf … cartesian components of force between tip and sample moderated by probe particle
4. OutX.xsf,OutY.xsf,OutZ.xsf … cartesian components of probe-particle position after relaxation
5. OutTSS1.xsf,OutTSS2.xsf … hopping from tip to probe-particle ( OutTSS1.xsf ) and form probe particle to atoms of surface ( OutTSS2.xsf, simple incoherent sum ). Both evaluated just as simple exponential
6. OutT1.xsf,OutT1.xsf … the same hoppings computed with angular dependence considerring two p-orbitals ( please see the paper supplementary or the code for more details )

1. Good software for visualization of .xsf files is VESTA
2. you can allso use my plotting scripts in python / numpy / pylab / matplotlib , however detailed documentation is beyon this tutorial. If you are familiar with python / numpy / pylab / matplotlib just look into STHM_Utils.py file where is most of important plotting rutines. The usage of plotting subroutines is obvious from plotAll.py script.
3. in order to obtain frequency shift you should convert OutFz.xsf file. In package of my python scripts it is done by subroutine 'deSared()' in STHM_Utils.py. For example in script plotAll.py is visible example of usage. There is an example of most basic stuff how to load vertiacal force, convert it to frequency shift and plot the images.

ilist  = range(30,50,2) # range of z-slice to plot ( make sure enought slices is in OutFz.xsf )
extent = (0,15,0,18)    # range of datafile in x,y direction
  
# imports
from STHM_Utils import *
from xsfutil import *
from pylab import *
   
# conversion from force to frequnecy shift
data_Fz,lvec, nDim, head = loadXSF('OutFz.xsf')   # this will load a .xsf file
df_data = F2df( data_Fz, dz = 0.1, k0 = 1800.0, f0=30300.0 ) # this will convert Fz to frequency shift with dz step of 0.1 Angsroem and cantilever stiffnes 1800 N/m and basic frequency 30.3 kHz
 
# plotting     
plotWithAtoms( df_data[:,1:-1,1:-1], ilist, extent, dz = 0.1,  atoms=None, cmap='gray' );     #                        
savefig('df.png',          bbox_inches='tight', pad_inches=0 )
show()