This tutorial explains how to perform scanning tunneling microscopy simulation based on the Tersoff-Hamann theory by using Al(111)(3x3)-O as an example. In this example, instead of using a simple local density of states from those obtained in a standard SCF calculations, the wave functions are reconstructed in such a way that they decay exponentially in the vacuum region (zero bias).
Suppose the structure is fully optimized and a new input file is created by using the geom2nfinp program in the utility directory.
0 0 0 0 0 0 : Al(111) (3x3) 5.0 15.0 2 37 37 : GMAX GMAXP KTYP KATM KATM2 1 0 : NUM_SPACE_GROUP TYPE Cartesian 16.235313117399 0.000000000000 0.000000000000 8.117656558700 14.060193598063 0.000000000000 0.000000000000 0.000000000000 44.186925502159 4 4 1 1 1 1 : K_MESH 1 0 : NCORD NINV : IWEI IMDTYP ITYP 2.705826513042 1.562208448277 1.542904250830 1 1 2 -0.101040174364 -0.058336519962 0.061753986538 1 1 1 5.512784204974 -0.058313698882 0.061749232330 1 1 1 10.823532872680 -0.010553333078 0.093326154025 1 1 1 2.705891678362 4.803367697679 0.061751269407 1 1 1 8.158736454981 4.710447914641 -0.150408894360 1 1 1 13.488359690248 4.710457156282 -0.150421161826 1 1 1 5.402626131859 9.378731841492 0.093325772329 1 1 1 10.823556250816 9.326033439020 -0.150420075404 1 1 1 16.244459949678 9.378743914505 0.093324680067 1 1 1 0.005685553104 3.121206464856 -4.442179390175 1 1 1 5.406084220392 3.121202931158 -4.442174570243 1 1 1 10.823537730940 3.082192093875 -4.504882024997 1 1 1 2.705880947950 7.818273168592 -4.436392497058 1 1 1 8.081025974571 7.832362395824 -4.504881761126 1 1 1 13.566048376540 7.832361823915 -4.504886108011 1 1 1 5.405657148113 12.494420562183 -4.436390486665 1 1 1 10.823544447840 12.504513279356 -4.442179755245 1 1 1 16.241420153507 12.494418723842 -4.436392516076 1 1 1 0.000000000000 -3.124487465505 -8.837385099216 1 0 1 5.411771037948 -3.124487465505 -8.837385099216 1 0 1 10.823542079420 -3.124487465505 -8.837385099216 1 0 1 2.705885520752 1.562243734036 -8.837385099216 1 0 1 8.117656558700 1.562243734036 -8.837385099216 1 0 1 13.529427596648 1.562243734036 -8.837385099216 1 0 1 5.411771037948 6.248974931245 -8.837385099216 1 0 1 10.823542079420 6.248974931245 -8.837385099216 1 0 1 16.235313117368 6.248974931245 -8.837385099216 1 0 1 0.000000000000 0.000000000000 -13.256077650288 1 0 1 5.411771037948 0.000000000000 -13.256077650288 1 0 1 10.823542079420 0.000000000000 -13.256077650288 1 0 1 2.705885520752 4.686731199776 -13.256077650288 1 0 1 8.117656558700 4.686731199776 -13.256077650288 1 0 1 13.529427596648 4.686731199776 -13.256077650288 1 0 1 5.411771037948 9.373462398956 -13.256077650288 1 0 1 10.823542079420 9.373462398956 -13.256077650288 1 0 1 16.235313117368 9.373462398956 -13.256077650288 1 0 1 13.0 0.5 26.982 1 1 0.0 : IATOMN ALFA AMION ILOC IVAN ZETA1 8.0 0.5 15.999 3 1 0.0 : IATOMN ALFA AMION ILOC IVAN ZETA1 0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC 0 1 : IPRE IPRI 200 1000 0 86400.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP 6 1 : MIXWAY MIXWHAT 0 50 0.30 : ITER_START KBXMIX MIX_ALPHA 0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM 200.00 2 1 0.10D-08 : DTION IMDALG IEXPL EDELTA -0.002 5.00D+03 0 : WIDTH FORCCR ISTRESS ggapbe 1 : XCTYPE KSPIN 1.00 : DESTM 102 : NBZTYP 4 4 4 : (DUMMY) 4 4 4 : (DUMMY) 68 : KEG 1 : NEXTST 0 : (DUMMY) 2 : IMSD 0 : EVALUATE_EKO_DIFF 0 : NPDOSAO 0 0.0 : (DUMMY) &ESM BOUNDARY_CONDITION BARE &ENDNote for the STM simulation (and electronic structure analysis such as density of states calculations), we recommend to use the Gamma-centered k-point mesh such as
4 4 1 1 1 1 : K_MESHNote also when we perform a STM simulation, make sure that ESM method (bare Coulomb (BARE) or BC1 for the boundary condition) is used, to obtain the accurate electrostatic potential in the vacuum region by setting.
&ESM BOUNDARY_CONDITION BARE &END
$ mpirun -np 24 < nfinp_scf > nfout_scf
Having obtained the self-consistent charge and potential, we then extract planar average of electrostatic potential by executing
$ state2chgpro.sh nfout_scf > chgpro.dat
To generate a gnuplot script to plot the planar average of charge density and potential, execute
$ chgpro2plt chgpro.dat > chgpro.gp
and type
$ gnuplot chgpro.gp
To perform the STM simulation, the vacuum level needs to be calculated. This can be done by first executing
$ chgpro2vac -stm
Then we obtained a file "stm" is generated, which looks like:
&STM STMOPT 1 Z0 14.11530000D0 VACUUM_LEVEL 0.00026382D0 &END
Here, Z0 indicate the z-coordinate at which the charge density is negligible small (default 1e-8 e/Bohr^3) and VACUUM_LEVEL is the vacuum level in Hartree. This is copy-pasted in the input file for the STM simulation as follows.
Having obtained self-consistent wave functions and vacuum level, we are ready to simulate the STM image.
0 0 0 0 0 0 : Al(111) (3x3) 5.0 15.0 2 37 37 : GMAX GMAXP KTYP KATM KATM2 1 0 : NUM_SPACE_GROUP TYPE Cartesian 16.235313117399 0.000000000000 0.000000000000 8.117656558700 14.060193598063 0.000000000000 0.000000000000 0.000000000000 44.186925502159 4 4 1 1 1 1 : K_MESH 1 0 : NCORD NINV : IWEI IMDTYP ITYP 2.705826513042 1.562208448277 1.542904250830 1 1 2 -0.101040174364 -0.058336519962 0.061753986538 1 1 1 5.512784204974 -0.058313698882 0.061749232330 1 1 1 10.823532872680 -0.010553333078 0.093326154025 1 1 1 2.705891678362 4.803367697679 0.061751269407 1 1 1 8.158736454981 4.710447914641 -0.150408894360 1 1 1 13.488359690248 4.710457156282 -0.150421161826 1 1 1 5.402626131859 9.378731841492 0.093325772329 1 1 1 10.823556250816 9.326033439020 -0.150420075404 1 1 1 16.244459949678 9.378743914505 0.093324680067 1 1 1 0.005685553104 3.121206464856 -4.442179390175 1 1 1 5.406084220392 3.121202931158 -4.442174570243 1 1 1 10.823537730940 3.082192093875 -4.504882024997 1 1 1 2.705880947950 7.818273168592 -4.436392497058 1 1 1 8.081025974571 7.832362395824 -4.504881761126 1 1 1 13.566048376540 7.832361823915 -4.504886108011 1 1 1 5.405657148113 12.494420562183 -4.436390486665 1 1 1 10.823544447840 12.504513279356 -4.442179755245 1 1 1 16.241420153507 12.494418723842 -4.436392516076 1 1 1 0.000000000000 -3.124487465505 -8.837385099216 1 0 1 5.411771037948 -3.124487465505 -8.837385099216 1 0 1 10.823542079420 -3.124487465505 -8.837385099216 1 0 1 2.705885520752 1.562243734036 -8.837385099216 1 0 1 8.117656558700 1.562243734036 -8.837385099216 1 0 1 13.529427596648 1.562243734036 -8.837385099216 1 0 1 5.411771037948 6.248974931245 -8.837385099216 1 0 1 10.823542079420 6.248974931245 -8.837385099216 1 0 1 16.235313117368 6.248974931245 -8.837385099216 1 0 1 0.000000000000 0.000000000000 -13.256077650288 1 0 1 5.411771037948 0.000000000000 -13.256077650288 1 0 1 10.823542079420 0.000000000000 -13.256077650288 1 0 1 2.705885520752 4.686731199776 -13.256077650288 1 0 1 8.117656558700 4.686731199776 -13.256077650288 1 0 1 13.529427596648 4.686731199776 -13.256077650288 1 0 1 5.411771037948 9.373462398956 -13.256077650288 1 0 1 10.823542079420 9.373462398956 -13.256077650288 1 0 1 16.235313117368 9.373462398956 -13.256077650288 1 0 1 13.0 0.5 26.982 1 1 0.0 : IATOMN ALFA AMION ILOC IVAN ZETA1 8.0 0.5 15.999 3 1 0.0 : IATOMN ALFA AMION ILOC IVAN ZETA1 21 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC 0 1 : IPRE IPRI 200 1000 0 86400.00 0 : NMD1 NMD2 LAST_ITER CPUMAX IFSTOP 6 1 : MIXWAY MIXWHAT 0 50 0.30 : ITER_START KBXMIX MIX_ALPHA 0.20 0.30 0.20 0.20 0.20 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM 200.00 2 1 0.10D-08 : DTION IMDALG IEXPL EDELTA -0.002 5.00D+03 0 : WIDTH FORCCR ISTRESS ggapbe 1 : XCTYPE KSPIN 0.10 : DESTM 102 : NBZTYP 4 4 4 : (DUMMY) 4 4 4 : (DUMMY) 68 : KEG 1 : NEXTST 0 : (DUMMY) 2 : IMSD 0 : EVALUATE_EKO_DIFF 0 : NPDOSAO 0 0.0 : (DUMMY) &ESM BOUNDARY_CONDITION BARE &END &STM STMOPT 1 Z0 14.11530000D0 VACUUM_LEVEL 0.00026382D0 BIAS 1.00 &ENDThis input file is exactly the same as that for SCF, except for
21 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACCwhere ICOND is set to 21 and the optional block
&STM STMOPT 1 Z0 14.11530000D0 VACUUM_LEVEL 0.00026382D0 BIAS 1.00 &ENDwhere Z0 and VACUUM_LEVEL are defined for the wave function reconstruction, as well as the STM bias in eV is set by using the keyword BIAS. The bias is also set by using the keyword DESTM in the main input file as (below XCTYPE)
0.10 : DESTM
$ mpirun -np 24 < nfinp_stm > nfout_stmand we obtain "nfchgt_r.data", which contains the local density of state (LDOS, or local charge) for the occupied (-DESTM - Fermi level to Fermi level) and unoccupied (from Fermi level to +DESTM - Fermi level) states, which is then used to obtained simulated images.
We are almost there -- Let us use a utility program "chg2xsf" to visualize the simulated data. First execute
$ chg2xsf -i nfinp_dav -c nfchgt_r.data -p stm
and we obtain LDOS files in the XSF format named "stm_occ.xsf" and "stm_unocc.xsf", respectively, which correspond to the constant height STM images.
The height can be changed in graphic softwares such as XCrySDen or VESTA.
Type
$ chg2xsf -h
for the usage of the chg2xsf program.
To obtain a simulated image corresponding to a constant current one, use
$ chg2xsf -i nfinp_dav -c nfchgt_r.data -stm 1.e-8 -p stm
Here by adding the option "-stm 1.e-8", the chg2xsf generates the simulated topographic image at the charge density of 1.e^-8 e/Bohr^3. Note the generated data (stm_occ.xsf and stm_unocc.xsf) are in 3D, but it is the repetition of 2D data in the z-direction (see the simulated image (to be added in this page)).
An STM simulation can be performed by using ICOND=10 (all the charge density) or ICOND=11 (soft part of the charge density). Try them and how the simulated images look like. Note by setting STMOPT=0 in the &STM...&END block, identical images with ICOND=11 should be obtained.