In this exercise we are going to see how we can generate Scanning Tunneling Microscope (STM) images from simulations. While this is a more elaborate example where we need additional scripts and techniques, the basic idea is simple: Imagine a STM cantilever scanning the surface in the constant current mode, the current is then proportional to the electron density, which is something we get from our simulation (remember the cube files?). Therefore, all we have to do is to find an isosurface of the electron density and plot this isosurfaces z-coordinate.

This exercise is based on ideas and material by Daniele Passerone.

• On the server is a package for you to unpack (hohoho ), containing a number of input files. Run the following in a new and empty directory:

  tar xf /users/tiziano/CHE437_ex7.tar.gz

• The scripts are contained in yet another python package:

  pip install --user https://github.com/ltalirz/asetk/archive/master.zip

  … and since you have setup the path variable in a previous exercise, you should now have the following new commands available: stm.py, cube-plot.py, cp2k-sumbias.py. If the installation fails, make sure that you do not have the CP2K module loaded: module list should return an empty list. To explicitly unload the CP2K module, run module unload cp2k.

First take a look at the complete geometry in all.xyz using vmd. As you can see, this is a much larger geometry than what we had before and would therefore take a lot longer to optimize. This is why some tricks are used in geo_opt.inp: For the nanoribbon (mol.xyz), we use the Density Functional based Tight-Binding (DFTB) method instead of GPW and couple the systems together by Force Fields (FIST). Look at the CP2K input file geo_opt.inp to get an idea how that looks like and then run the simulation.

Like with the geometry optimization in the previous examples, you now should have a file nanoribbon-geo_opt-pos-1.xyz. From this we have to extract only the molecule into a new xyz file since we are going to create a STM image only for the molecule, not the substrate:

echo "224" > nanoribbon.xyz  # the molecule contains 224 atoms, first line in a XYZ file
echo "optimized using DFTB/FF" >> nanoribbon.xyz  # second line in a XYZ file, an optional comment
tail -n 1064 nanoribbon-geo_opt-pos-1.xyz | head -n 224 >> nanoribbon.xyz  # the complete system contains 1064, but we want only the first 224 (= the molecule)

To get the actual electron density, we are now going to run a full DFT calculation using a large basis set (TZV2P) on the nanoribbon geometry. This time, the input nanoribbon.inp is similar to what you are used to, except for the fact that to speedup the calculation, we use the wavefunctions from the previous calculation as a starting point (the RESTART_FILE_NAME option).

This calculation will take a while to finish: run it in parallel using 4 processes (mpirun -np 4 …) and in the background.

At the end of it, you can visualize the orbitals using either VMD as shown in previous exercises, or using one of the scripts provided in the newly installed Python package, which will generate a series of images (nanoribbon-WFN_*.png), visualizing the different orbitals:

cube-plot.py --cubes *WFN*.cube --plot --position 22  # the position for the contour plot should be 2 angstrom above the surface

To get an actual STM image, we now have to combine the wavefunctions into a single one:

# use your output file of the full DFT calculation as your levelsfile!
cp2k-sumbias.py --cubes *WFN*.cube --levelsfile nanoribbon.out --vmin -2.0 --vmax 2.0 --vstep 0.5 | tee sumbias.out
# and pipe the output to the file sumbias.out and the screen simultaneously by using 'tee'

The parameters –vmin, –vmax and –vstep determine which bias voltages for the tip (the potential between the substrate/molecule and the tip) you want to simulate, in our case $-2.0$, $-1.5$, … $2.0$.

It is important to note that for a given bias voltage, for example $-2.0$ (current goes from the substrate/molecule to the tip) all orbitals with an energy between $-2.0 eV$ and $0 eV$ have to be taken into account.

At this point you should have a new set of combined CUBE files: stm_-2.00V.cube..stm_+0.00V.cube..stm_+2.00V.cube, one for each bias voltage, containing the respective electron density.

From these we can finally generate the actual STM images, which should give you a set of files stm_*V.cube.iso1e-07.png:

# zcut is the minimum z-height
stm.py --stmcubes stm_*.cube --isovalues 1.0e-7 --zcut 22 --plot

Why are there no images for certain bias voltages? Would you expect the same for a metallic substrate?