:sequential_nav: next .. _tutorial-stm-images: Simulation of STM images ============================ :Author: Jorge Pilo (ICN2) .. sidebar:: **Have you set up the local environment?** If not, :ref:`do that now ` before proceeding. This tutorial provides some guided examples for the use of the tools available to simulate STM and STS images with Siesta. The important tools are two: the ``plstm`` executable and the ``stm`` executable (the latter one is also called ``ol-stm`` in this tutorial and it is also known with the more appropriate name `wfs2ldos`). Before starting it is important to read :ref:`here` to get an overview of the tools. .. note:: For some background and a gallery of images, see the `presentation `_ by Jorge Pilo. STM image of a benzene molecule ------------------------------- In this exercise we will obtain the STM image of a benzene molecule, We will compare the ``plstm`` and ``ol-stm`` ways and see how these tools can complement each other. Options in the Siesta run ......................... First of all, enter the folder ``STM_Benzene`` and explore the `benzene.fdf` file. The keywords we should focus on are: * The **%block LocalDensityOfStates** that asks Siesta to calculate the LDOS in an energy range. * **write-denchar** that asks Siesta to generate two files called `SystemLabel.PLD` and `SystemLabel.DIM` where information needed to plot the charge density and/or wavefunctions in real space is dumped. * **COOP.Write** the writes information on the wavefunctions at all k-points used to sample the BZ during self-consistency. * **SaveElectrostaticPotential** to save the Hartree potential * **WFS.Energy.Min** and **WFS.Energy.Max** to limit the number of wave functions plotted to the ones corresponding to an eigenvalue on this energy range. .. note:: Some considerations about the structure: we can study only planes perpendicular to the z direction, therefore the structure must be prepared accordingly. Run the calculation with: .. code-block:: bash siesta < benzene.fdf > benzene.out STM image from direct LDOS calculation: plstm ............................................. We first use ``psltm`` to obtain a first STM image. For this tool only the **%block LocalDensityOfStates** was necessary in the siesta input and only the `benzene.LDOS` file is needed. This file contains LDOS(x,y,z) in the same grid used by Siesta in that calculation (i.e., it will depend on the MeshCutoff used). This LDOS function is periodic. Create a new folder and copy or link there the relevant file: .. code-block:: bash mkdir PLSTM cd PLSTM ln -sf ../benzene.LDOS ./benzene.LDOS To get the STM image at a distance of 9.82 Bohr from z=0: .. code-block:: bash plstm -z 9.82 -o stm_9.82bohr benzene.LDOS The "stm9.82bohr" is the output file that contains the STM data plottable with gnuplot. Open gnuplot (just type ``gnuplot``) and enter the commands: .. code-block:: gnuplot set pm3d set pm3d map set xrange [0:22.253958] set yrange [0:20.845341] set palette rgb 34,35,36 set terminal wxt size 440,440 splot 'stm_9.82bohr' Congratulations! You got your first STM image! To exit gnuplot type ``quit``. More on the ``plstm`` will be cover in the next sections of this tutorial. STM image from wavefunction information using ol-stm .................................................... The ``plstm`` program is very simple to use but it is limited in flexibility. For instance it does not allow to change the energy range or the grid for the calculation (unless you run again a new Siesta calculation). A more flaxible tool is the ``ol-stm`` program (called with the ``stm`` executable for historic reasons). The ``ol-stm`` starts from the wavefunctions to first generate the LDOS, from where the STM/STS images can be obtained. Importantly, this code can optionally compute the values of the wavefunctions generated by Siesta on a reference plane and then extrapolate the value of these wavefunctions into the vacuum. In fact, for realistic experimental conditions the STM tip is relatively far from the surface. At this distance the accuracy of the Siesta calculation is not sufficiently good, because the basis set is short-ranged (and also because of other issues common to DFT approaches, including plane-wave simulations). To be able to extrapolate the value of the wavefunctions is a fairly important feature. More can be understood reading :ref:`here`. Let's go back to the folder of the siesta calculation ``cd ..``. Located in the folder one can see the files generated by the Siesta run: * `benzene.PLD` * `benzene.DIM` * `.ion` * `benzene.fullBZ.WFSX` * `benzene.VH` All these files are needed by the ``stm`` executable. .. warning:: ``stm`` imposes to have the wavefunctions written in the `SystemLabel.WFSX`. Siesta returns this info in different files depending on the way used to obtain them. In our case we have the `benzene.fullBZ.WFSX` file name. This must be changed to `benzene.WFSX`. .. code-block:: ln -sf benzene.fullBZ.WFSX benzene.WFSX The file `benzene-stm.fdf` was instead already present in the folder and it is the input for the ``stm`` executable. Explore this file. A part from the system information, some important keywords are: * **STM.MinZ** (real length): Defines the minimum value of the z-component of the 3D grid. * **STM.MaxZ** (real length): Defines the maximum value of the z-component of the 3D grid. * **STM.RefZ** (real length): Defines the z-component position of the reference plane from where the wavefunctions will be extrapolated into vacuum. The LDOS values for planes below the reference level are obtained from the original wavefunctions. Those above the reference level are obtained from projected wavefunction values. * **STM.VacZ** (real length): Defines the z-component position of the plane on which the zero of the Hartree potential in the vacuum region is taken. This should be a position at the center of the vacuum of the slab. You might want to check (better: plot) the V(z) values printed in the `benzene.v_ave_z` file. This file is produced by ``stm`` itself. Therefore you need a first guess and maybe rerun the calculation in case this value was wrong. * **STM.Emin** (real energy): Lower value of the energy window (eV). * **STM.Emax** (real energy): Upper value of the energy window (eV). * **STM.NumberPointsX** (integer): Number of subdivision of the grid in the direction of the first lattice vector. * **STM.NumberPointsY** (integer): Number of subdivision of the grid in the direction of the second lattice vector * **STM.NumberPointsZ** (integer): Number of subdivision of the grid in the z-direction (which must be the direction of the third lattice vector) More keywords can be explored in the :ref:`reference manual` for the `stm` (or `wfs2ldos` program). We can now run .. code-block:: bash stm < benzene-stm.fdf > benzene-stm.out The file generated are: * `benzene.CELL_STRUCT`: Cell structure parameters * `benzene.STM.LDOS`: In STM mode, it contains the value of the LDOS (all relevant spin components) at the grid defined by the user, in Siesta grid format. * `benzene.v_ave_z`: An ASCII file containing the V(z) profile of the electrostatic potential. Plotting the information in this file is useful to set appropriately STM.Vacz and STM.RefZ. To obtain the STM data we are going to use the ``rho2xsf`` util. We need only the `benzene.XV` file (produced by any Siesta run and containing the structure) and the file with the LDOS. Better create a separate folder, and take into account that the ``rho2xsf`` recognizes only files with single extension, therefore: .. code-block:: bash mkdir STM cd STM ln -sf ../benzene.XV ./benzene.XV ln -sf ../benzene.STM.LDOS ./benzene.LDOS Type ``rho2xsf``. The program opens an interactive shell and you have to insert the following info: [Maybe better to prepare an input file to represent the standard input stream] .. code-block:: Specify SystemLabel (or ’siesta’ if none): benzene Would you use Bohr (B) or Ang (A) ? B Enter reference point in Bohr : 0.0 0.0 0.0 Enter 1st spanning vector in Bohr : 22.253958475 0.000000000 0.000000000 measured from reference point in Bohr : 0.000000000 22.253958475 Enter 2nd spanning vector in Bohr : 0.000000000 20.845341211 0.000000000 measured from reference point in Bohr : 0.000000000 20.845341211 Enter 3rd spanning vector in Bohr : 0.000000000 0.000000000 37.794537555 measured from reference point in Bohr : 0.000000000 37.794537555 Enter number of grid points along three vectors 40 40 40 Add grid property (LDOS, RHO, ...; or BYE if none): LDOS Add grid property (LDOS, RHO, ...; or BYE if none): BYE The benzene.XSF file has been created and can be opened with the ``xcrysden`` program. Type ``xcrysden``. Go to **File** -> **Open Structure** -> **Open XSF**, select the XSF file just generated and "OK". Then go to **Tool** -> **Data Grid**. The program should already select the LDOS file so just "OK". A new window is opened to select the isosurface. Insert 0.01 as the value of the isosurface and untick the "Display isosurface" on the top left of the window. Go to the "Plain #1" tab. "Select color basis" to "BLACK-BROWN-WHITE", tick "display color-plane" and "Submit". Pictures of the process are in the presentation mentioned at the beginning of this tutorial. You have the STM of the benzene. .. note:: An alternative would be to create a .cube file where the STM image is represented. This is done using ``g2c_ng``. .. code-block:: bash g2c_ng -s benzene.CELL_STRUCT -g benzene.STM.LDOS And the .cube file can be visualized using VESTA (not by default in the Siesta Mobile) A variety of images of the benzene STM can be found in the slide presentation mentioned at the top of this tutorial. .. note:: One can now use ``plstm`` to postprocess the `benzene.STM.LDOS` file and get the appropriate slices (either constant-current or constant-height). The results can be plotted with gnuplot. STM images for a monolayer of Cr with Spin-Orbit coupling --------------------------------------------------------- We will now explore the possibility to plot STM images for magnetic systems. Enter the folder STM_Cr_3atms You can see that the `Cr.fdf` file contains the keywords for the initialization of the spin and to activate spin orbit coupling (if you are not familiar with these concepts look at the Spin-orbit coupling tutorial). Run siesta, it will take about 5 minutes: .. code-block:: bash siesta < Cr.fdf > Cr.out Is important to check that you are in good agree with the final magnetic moments direction, these are shown because we specified the keyword **WriteMullikenPop 1** in the input file. We show you the final directions for this example (can be found in the "mulliken" section of the `Cr.out`): .. code-block:: bash ATOM Charge Mag.Mom. Sx Sy Sz 1 5.99983 4.93634 2.466 4.276 0.000 2 6.00020 4.93714 2.470 -4.275 -0.000 3 5.99997 4.93661 -4.937 0.003 -0.000 Total 18.00000 0.00406 -0.001 0.004 -0.000 It is important to say that these values are in cartesian coordinates. Get the stm image with ``plstm``: .. code-block:: bash mkdir PLSTM cd PLSTM ln -sf ../Cr.RHO ./Cr.RHO plstm -z 4.46 -s s -X 2 -Y 2 -o stm_s_4.46bohr Cr.RHO To notice: * We now used the .RHO file, in fact the ``plstm`` recognizes both LDOS and RHO, so they can be used interchangeably. * We used the spin option ``-s`` set to ``s`` to plot the total spin magnitude [Spin(up)-Spin(down)] instead of the total charge (``q`` that is the default). * We used -X and -Y to ask to plot more than just one unit cell. You can use the `School_STM.gplot` script to plot the result with gnuplot: .. code-block:: bash ln -sf ../School_STM.gplot ./School_STM.gplot gnuplot School_STM.gplot Run with no spin: .. code-block:: bash plstm -z 4.46 -s q -X 2 -Y 2 -o stm_q_4.46bohr Cr.RHO and selecting just one spin component: .. code-block:: bash plstm -z 4.46 -s x -X 2 -Y 2 -o stm_x_4.46bohr Cr.RHO Remember to change the file name inside the gnoplot script in order to plot the new files. Images with magnetic tips ......................... The use of the spin option of the ``plstm`` (``-s``) is useful to explore the spin dependence of the LDOS, however a scenario that is more close to an experimental set-up is the situation when a magnetic tip is used. We propose two different vectors for the magnetic moment of the tip: V1=(0.5, 0.866, 0) and V2=(0.866, -0.5, 0). Important to say that these vectors have to be normalized. Type: .. code-block:: bash plstm -z 4.46 -v '0.5 0.866 0' -X 3 -Y 3 -o stmV14.46bohr Cr.RHO And then: .. code-block:: bash plstm -z 4.46 -v '0.866 -0.5 0' -X 3 -Y 3 -o stmV24.46bohr Cr.RHO Use the script `School_STM.gplot` to post processing both files, after changing the file name inside the script. Look at the presentation at the beginning of this page for more insights on the results. If we are not happy with the grid or energy we can use the procedure involving the ``ol-stm`` program explained in the previous section. More on that on the slides at the beginning of this page. Conclusions ----------- We presented some simple examples to show how to use the STM utils from Siesta. In a real life calculation it is important to be sure about the convergence of all the parameters and inputs that the job needs. As we could see, there are different ways to porecess a STM image with SIESTA, and this overview is not sufficient to cover all the aspects of the topic. For more info you can check the user’s manual for the codes (links from :ref:`here`). There are some articles in which the authors use the STM Tool from SIESTA. We invite you to read: * \E. Cisternas and J.D. Correa: *Theoretical reproduction of superstructures revealed by STM on bilayer graphene*, Chemical Physics 409, 74-78 (2012). DOI: `10.1016/j.chemphys.2012.09.021 `_. * Ana Sanz-Matías, *et al.*: *Computational insight into the origin of unexpected contrast in chiral markers as revealed by STM*, Nanoscale 10, 1680-1694 (2018). DOI: `10.1039/C7NR07395J `_.