Creating AFM and STM images

This is a short workflow regarding the creation of STM and AFM images with different quantities outputted in .cube files from the FHI-aims code. Scripts and older documents that aid visualization are also linked here.

STM simulations

FHI-aims and many other electronic structure codes can easily produce the necessary ingredients to simulate STM images within the Tersoff-Hamann approximation. Within the framework of DFT and under some applied voltage bias \(V\), the tunneling current shows the following proportionality

\[\begin{split}I_t(\boldsymbol{r}) \propto \begin{cases} \sum_{\boldsymbol{k}, \nu} |\psi_{\boldsymbol{k},\nu}(\boldsymbol{r})|^2 f(\varepsilon_{\boldsymbol{k}, \nu} - E_f) [1-f(\varepsilon_{\boldsymbol{k}, \nu} - (E_f-eV))] & \text{if } V < 0 \\ \sum_{\boldsymbol{k}, \nu} |\psi_{\boldsymbol{k},\nu}(\boldsymbol{r})|^2 [1-f(\varepsilon_{\boldsymbol{k}, \nu} - E_f)] f(\varepsilon_{\boldsymbol{k}, \nu} - (E_f-eV)) & \text{if } V > 0 \end{cases}\end{split}\]

where \(\psi_{\boldsymbol{k}, \nu}\) are the Kohn-Sham orbitals with respective energies \(\varepsilon_{\boldsymbol{k}, \nu}\) and \(E_f\) is the Fermi energy. The condition on the occupation functions make sure that the sums only include states within an energy window between \(E_f\) and \(V\). Note that if \(V>0\), the current probes unoccupied states. \(\sum_{\boldsymbol{k}, \nu} |\psi_{\boldsymbol{k},\nu}(\boldsymbol{r})|^2 \delta(\epsilon - \varepsilon_{\boldsymbol{k}, \nu})\) is usually called the local density of states LDOS.

As explained in the manual of FHI-aims, once one has a system geometry for which an STM image should be created, a full SCF convergence including the control.in command needs to be run:

output cube stm 0.3

Note that the last value is given in eV and it can be negative or positive. It corresponds to the energy window, starting from the Fermi energy determined in the calculation. At the end of this calculation, a cube file will be generated where the LDOS is given at each space point in the grid.

There are two possible modes of producing STM images: constant current and constant height.

For the constant height image, different values of current will be present at each position, and building an STM image corresponds to plotting all values of the cube file that lie at a given distance \(d\) above the sample of interest. The color code is the value of the LDOS at each point with a fixed d. These images are easily created by just taking a slice at any given height in the cube file that was generated.

For the constant current image, building an STM image corresponds to finding an isosurface in the cube file that maintains the value of the LDOS constant, within a certain threshold. The color code is the height of each point for a fixed value of the LDOS.

Because building an image in the constant current case is a bit more involved, but also more common to be of interest, we provide a script and a few steps to make such images below.

The idea is that, once the cube file is generated:

  • Find coordinates of each point and its corresponding STM current

  • Interpolate those points

  • Find the coordinates of a current isosurface

  • Plot these points on a plane (normally xy) colored by height (normally the z coordinate)

Scripts that do these procedures can be found here. For historical reasons, there are two scripts, but making them a single workflow is in the makings (should someone volunteer).

The first script: position-voltage-extract.py will extract the Cartesian coordinates and the current of each point and write a file with four columns.

The second script: plot_stm_images.py will then do all the following steps, taking the file generated by the first script as an input. The help function is here for completeness:

plot_stm_images.py -h
usage: plot_stm_images.py [-h] [--iso-value ISO_VALUE] [--z-min Z_MIN] [--z-max Z_MAX] [--scatter-output SCATTER_OUTPUT] [--density-output DENSITY_OUTPUT] input_file

Interactive STM data analysis script

positional arguments:
   input_file            input data file

options:
  -h, --help            show this help message and exit
  --iso-value ISO_VALUE
                    isosurface value
  --z-min Z_MIN         Minimum z-coordinate
  --z-max Z_MAX         Maximum z-coordinate
  --scatter-output SCATTER_OUTPUT
                    output file for scatter plot
  --density-output DENSITY_OUTPUT
                    output file for density plot

Feel free to play with the images, isovalues, etc.

AFM simulations

For the AFM simulations, we rely on the probe-particle code developed by the Nanosurf Lab in Prague.

In this case, once one has a system geometry for which an AFM image should be created, a full SCF convergence including the control.in command needs to be run:

output cube hartree_potential

For FHI-aims data, for molecules on surfaces, we have experienced a better performance when we removed the surface and the lattice periodicity (ensuring the molecules are not cut through the periodic boundaries) before generating these files. Of course, in principle it should work also with periodicity, but the parameter range seems difficult to optimize. We will update this page with solutions when we find them.

If running the “non-periodic” case, an example of the input file for the ppafm code v.0.4.0 (for an image without any V bias and adapted from their website) and the commands that need to be run follow:

params.ini:

probeType      54                      # atom type of ProbeParticle (to choose L-J potential ),e.g. 8 for CO, 54 for Xe
charge         +0.3                     # effective charge of probe particle [e]
klat            0.2                             # [N/m] harmonic spring potential (x,y) components, x,y is bending stiffnes
krad            20.00                           # [N/m] harmonic spring potential R component, R is the particle-tip bond-length stiffnes
r0Probe         0.0 0.0  4.00                   # [Å] equilibirum position of probe particle (x,y,R) components, R is bond length, x,y introduce tip asymmetry
PBC             False                           # Periodic boundary conditions ? [ True/False ]
gridN           642 213 874             # grid points for x, y, z (scaled from lattice dimensions)
gridA           47.3   0.0 0.0  # a-vector (x-dimension lattice vector)
gridB           0.0 28.6 0.0  # b-vector (y-dimension lattice vector)
gridC           0.0 0.0 41.7  # c-vector (z-dimension lattice vector)
scanMin         -3.6 -3.6       47            # start of scanning (x,y,z)
scanMax         43.5  23.3      58          # end of scanning (x,y,z)
scanStep        0.10 0.10 0.10
Amplitude       2.0             # [Å] oscilation amplitude for conversion Fz->df

commands:

ppafm-generate-elff -i cube_001_hartree_potential.cube
ppafm-generate-ljff -i cube_001_hartree_potential.cube
ppafm-relaxed-scan --Vbias 0.0 --pol_t 1.0 --pol_s 1.0
ppafm-plot-results --Vbias 0.0 --save_df --df

TODO: fix links to scripts and add actual examples of AFM and STM images