====== Simulation of STM images for a graphene nanoribbon adsorbed on a metallic substrate ======

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 [[exercises:2017_ethz_mmm:stm|ideas and material by Daniele Passerone]].


===== Preparation =====

  * 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: <code>tar xf /users/tiziano/CHE437_ex7.tar.gz</code>
  * The scripts are contained in yet another python package: <code>pip install --user https://github.com/ltalirz/asetk/archive/master.zip</code>... and since you have setup the path variable in [[exercises:2017_uzh_cmest:phonon_calculation|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''.

===== Geometry optimization =====

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:

<code bash>
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)
</code>

===== Calculating the nanoribbon =====

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:

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

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

<code bash>
# 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'
</code>

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'':

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

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