Part 3: STM simulations¶
By the end of this tutorial, you will be able to:
- Calculate a partial charge density of a graphite slab using VASP postprocessing.
- Simulate and plot a simulated constant-height scanning-tunneling-microscopy (STM) image at several tip heights.
9.1 Task¶
Perform a static calculation of a graphite slab to converge the charge density. Use the resulting WAVECAR file as input for a post-processing step to compute the charge density of bands just below the Fermi energy and plot a simulated STM image.
A scanning-tunneling-microscope (STM) utilizes quantum tunneling to image a clean and (semi)conducting surface on the atomic scale. One can measure the current at a constant tip height or track the tip position while keeping the tunneling current constant. The measured data is a function of the tip position, the bias voltage, and the local density of states. Thus, it is not always obvious what surface feature is represented by bright or dark regions in the image. Simulating STM images is an effective way to interpret experimental results and gain more insight into surface structure, and reconstructions [J. Phys. Chem. C, 2021, 125, 48, 26711–26717], and even dynamics [Nat. Mater, 2016, 15, 450–455].
A simple and effective way to simulate STM pictures with DFT calculations is the Tersoff-Hamann approach [Phys. Rev. Lett., 1983, 50, 1998], which is based on Bardeen's tunneling theory [Phys. Rev. Lett., 1961, 6, 57]. It models the tip as a structureless point corresponding to a single s-like wavefunction. This untangles the contributions to the tunnel current of the surface from the contributions of the tip, which are neglected. The approximation is valid for many surfaces, especially at larger distances, and is implemented in py4vasp. More involved approaches, like Chen's derivative rules [Phys. Rev. B, 1990, 42, 8841], are sometimes necessary for accurate results [J. Phys. Chem. C, 2021, 125, 48, 26711–26717], but require explicit treatment of the tip.
This exercise simulates an STM image of graphite in constant height mode with a bias voltage of -250 mV. This is achieved by plotting the smoothed partial charge density at a fixed distance above the surface.
9.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e09_STM_of_graphite. Check them out!
Graphite slab
2.46
1.0000000000000000 0.0000000000000000 0.0000000000000000
-0.5000000000000000 0.8660254037844386 0.0000000000000000
0.0000000000000000 0.0000000000000000 10.9105691056910570
C
10
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.3333333333333333 0.6666666666666666 0.0000000000000000
0.3333333333333333 0.6666666666666666 0.1250000000000000
0.6666666666666666 0.3333333333333333 0.1250000000000000
0.0000000000000000 0.0000000000000000 0.2500000000000000
0.3333333333333333 0.6666666666666666 0.2500000000000000
0.3333333333333333 0.6666666666666666 0.3750000000000000
0.6666666666666666 0.3333333333333333 0.3750000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.3333333333333333 0.6666666666666666 0.5000000000000000
SYSTEM = Graphite slab
ENCUT = 400
ISMEAR = 0
SIGMA = 0.2
ALGO = Fast
EDIFF = 1.0E-6
NEDOS = 3001
LORBIT = 11
SYSTEM = Graphite slab partial charges
ENCUT = 400
LPARD = .TRUE.
LPARDH5 = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3
EINT = -0.25 0.0
K-Points
0
Gamma
17 17 1
0 0 0
Pseudopotential of C.
The POSCAR file describes 5 layers of graphite using the experimental lattice parameter and interlayer spacing.
The INCAR file in $TUTORIALS/surface/e09_STM_of_graphite creates a simple, non-spin-polarized, static calculation, with increased DOS evaluation points and tightened convergence criterion.
Generally, it is recommended to set PREC = Single for STM simulations. This ensures that there is no interpolation of the charge on the fine FFT grid, since NGX = NGXF (equivalent for NGY and NGZ). The interpolation can result in high-frequency oscillations of the charge density in the vacuum region. py4vasp uses Gaussian smoothing to remove those potential oscillations when simulating STM pictures, so we keep the default PREC = Normal, doubling the resolution of our STM picture NGX times NGY to NGXF times NGYF.
The KPOINTS and POTCAR files are also standard and identical in both the main folder of exercise e09 and the subfolder for the partial charge density.
The second INCAR in $TUTORIALS/surface/e09_STM_of_graphite/partial_charges is used for the post-processing calculation of the partial charge density. Please read the section band-decomposed charge densities in the wiki and also the entries for the given INCAR tags LPARD, LPARDH5, LSEPK, LSEPB, NBMOD, and EINT to understand them. Make comments after each tag in the INCAR file reflecting what you learned.
Click to see the answer!
The LPARD tag toggles the partial charge postprocessing on or off.
LPARDH5 redirects the output from one or several PARCHG files to the vaspout.h5 hdf5 file.
LSEPK and LSEPB are both set to .FALSE. to sum up the output for all considered k points and bands instead of separating them.
EINT defines the energy interval for which the partial charge density is written, and NBMOD = -3 ensures that this interval is interpreted as relative to the Fermi energy.
These settings select the charge density of all bands and at all k points with eigenvalues between $\epsilon_F - 0.25$ eV and the Fermi energy $\epsilon_F$. This corresponds to a negative bias voltage of ~0.25V.
9.3 Calculation¶
Open a terminal, navigate to this example's directory, and run the self-consistent VASP calculation:
cd $TUTORIALS/surface/e09_*
mpirun -np 2 vasp_std
Use py4vasp to visualize the DOS close to the Fermi energy. This is the region where we will extract the partial charge density needed for the STM simulation.
Would the results differ when a small positive bias voltage is used instead of a small negative one?
import py4vasp
calc = py4vasp.Calculation.from_path( "./e09_STM_of_graphite/")
calc.dos.plot(selection="9,10")
Click to see the answer!
Since the DOS is symmetric in the region of $\sim \epsilon_F \pm 0.5$ eV, we expect nearly identical images for small positive or negative bias voltages. However, since the total DOS cannot tell us anything about the local real-space distribution of the contributing orbitals, we should look at the $p_z$ orbital contributions, which should contribute the most to the charge density normal to the surface. They are also symmetric around $\epsilon_f$ in the relevant energy range, as can be seen by adding the selection="pz" keyword argument to the plot method of the dos object.
Generally, a negative bias voltage measures the occupied states, while a positive voltage probes the unoccupied bands.
Copy the WAVECAR file to the partial_charges subdirectory, navigate there, and rerun VASP. (This is a quick post-processing step and is not parallelized. There is no need to run with mpi, although it will work.)
Take note of the output on your terminal. Which k points and bands contribute to the narrow energy window we selected?
cp $TUTORIALS/surface/e09_*/WAVECAR $TUTORIALS/surface/e09_*/partial_charges/.
cd $TUTORIALS/surface/e09_*/partial_charges/
vasp_std
Click to see the answer!
Only 2 bands, numbers 18 and 19 contribute in total, and both only contribute at k point number 33.
Use py4vasp to plot the simulated STM picture. Compare it with the STM image of graphite below. It was measured at ambient conditions (at room temperature, and in air), at the Department for Earth and Environmental Sciences, LMU and Center for NanoScience (CeNS), Munich. Images measured at low temperatures and in a vacuum are usually much clearer and less noisy. But even here we see atomic resolution and a clear triangular grid.
import py4vasp
calc = py4vasp.Calculation.from_path( "./e09_STM_of_graphite/partial_charges")
calc.partial_density.to_stm(supercell=7, tip_height=3)
The simulated image produces a triangular grid and matches the real image well if accounting for an arbitrary rotation. The surface atom with the direct coordinates 0.0 0.0 0.5 is visible as a bright spot, while the other surface atom at 0.333 0.666 0.5 only produces a small signal. Feel free to modify the tip height and the bias voltage (using the EINT tag in the post-processing step) to see how the simulated image changes.
If you use a larger negative bias voltage, more bands will contribute. The second surface atom of graphite, which has another carbon atom beneath it in the subsurface layer, will then become visible. The triangular lattice will be transformed into a hexagonal one.
9.4 Questions¶
- Can we simulate the STM tip properties with py4vasp's STM simulation implementation?
- Why is NBMOD = -3 a convenient setting for STM simulations and how does it relate to the bias voltage in the experiment?
- Why does the second surface atom only show up in the simulated image if the magnitude of the bias voltage is increased?
By the end of this tutorial, you will be able to:
- Calculate a partial charge density of a graphene sheet using VASP postprocessing.
- Simulate and plot a constant-current scanning-tunneling-microscopy (STM) image.
10.1 Task¶
Perform a static calculation of a graphene sheet to converge the charge density. Use the resulting WAVECAR file as input for a post-processing step to compute the charge density in a small interval around the Fermi energy and plot a simulated constant-current STM image.
A scanning-tunneling-microscope (STM) utilizes quantum tunneling to image a clean and (semi)conducting surface on the atomic scale. More information on STM methodology can be found in Example 9.
In constant current mode, the height of the STM tip is adjusted during the lateral scan so that the tunneling current remains constant. The height at each lateral point then forms the image. This is the most common mode of STM data acquisition.
In py4vasp constant current mode is turned on by setting selection="constant_current" in the to_stm method of the partial_charge class, while the default is the "constant_height" mode. A tunneling current in nA is then selectable. It defaults to 1 nA. The simulated picture is created by scanning down toward the surface from the middle of the vacuum for each lateral grid point until a value of the smoothed partial charge density, proportional to the selected current, is reached.
10.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e10_STM_of_graphene. Check them out!
Graphene sheet
2.46
1.0000000000000000 0.0000000000000000 0.0000000000000000
-0.5000000000000000 0.8660254037844386 0.0000000000000000
0.0000000000000000 0.0000000000000000 4.8780487804878050
C
2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.3333333333333333 0.6666666666666666 0.0000000000000000
SYSTEM = Graphene sheet
ENCUT = 400
ISMEAR = 0
SIGMA = 0.2
ALGO = Fast
EDIFF = 1.0E-6
NEDOS = 3001
LORBIT = 11
SYSTEM = Graphene sheet partial charges
ENCUT = 400
LPARD = .TRUE.
LPARDH5 = .TRUE.
NBMOD = -3
EINT = -0.25 0.0
K-Points
0
Gamma
17 17 1
0 0 0
Pseudopotential of C.
The POSCAR file only has one layer of carbon atoms, but the same lattice parameters as graphite in Example 9.
The KPOINTS and INCAR files are also identical to Example 09. Remind yourself of the meaning of all tags in the INCAR file and add comments describing their purpose.
Why are the tags LSEPK and LSEPB not set in the $TUTORIALS/surface/e10_STM_of_graphene?partial_charges/INCAR file?
Click to see the answer!
The default values of both LSEPK and LSEPB are .FALSE.. Thus we do not necessarily have to specify them. It is nevertheless a good idea to do so since it avoids ambiguity.
py4vasp will only allow you to plot an stm image simulation if neither LSEPK nor LSEPB are set to .TRUE.
10.3 Calculations¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e10_*
mpirun -np 2 vasp_std
This will run a static minimization of the electronic structure and prepare a WAVECAR file for post-processing. It will also allow you to visualize the DOS of graphene with py4vasp and compare it to that of the top layer of the graphite slab from the last exercise. Note that you can de-select the total DOS of both calculations by clicking on the labels on the right of the plot. Zoom into the region around the Fermi energy. The famous zero-gap electronic structure of graphene is not reproduced correctly due to the k point mesh not being dense enough in this example.
How do the two DOSs differ close to the Fermi energy?
import py4vasp
calc_graphene = py4vasp.Calculation.from_path("./e10_STM_of_graphene/")
calc_graphite = py4vasp.Calculation.from_path("./e09_STM_of_graphite/")
#calc_graphene.dos.plot(selection="1,2").label("graphene") + calc_graphite.dos.plot(selection="9,10").label("graphite")
a = calc_graphene.dos.plot(selection="1,2").label("graphene").to_plotly()
b = calc_graphite.dos.plot(selection="9,10").label("graphite").to_plotly()
# Add all traces from b into a
b.data[0]['line'].color = "#A82C35"
b.data[1]['line'].color = "#89ad01"
b.data[2]['line'].color = "#424242"
for trace in b.data:
a.add_trace(trace)
# Optionally update layout: add title, legend, etc.
a.update_layout(title="DOS Graphene and Graphite",
legend_title_text="Selection")
a.show()
Click to see the answer!
In the graphene sheet, both atoms have the equivalent DOS, which lies between the graphite DOS of atom 9 and atom 10.
The famous zero-gap electronic structure of graphene is not quite reproduced in this calculation. It would require significantly more k points to close the finite gap of about 1 eV. For STM picture simulation, however, this is of little significance.
Now let us compute the partial charge close to the Fermi energy of graphene ($\epsilon_F \pm 0.1$ eV).
Copy the WAVECAR file to the partial_charges subdirectory, navigate there, and run VASP again.
cp $TUTORIALS/surface/e10_*/WAVECAR $TUTORIALS/surface/e10_*/partial_charges/.
cd $TUTORIALS/surface/e10_*/partial_charges/
mpirun -np 2 vasp_std
10.5 Simulating and plotting the STM image¶
Use py4vasp to simulate and plot a constant current STM image. The negative bias voltage is again 250meV, to facilitate a good comparison with Example 9. A small bias voltage is appropriate for freestanding graphene, which would otherwise experience significant attraction to the tip [J. Vac. Sci. Technol. B, 2013, 31, 04D103].
py4vasp supports the selection of spin channels ("up", "down", or "total") which could lead to different images for magnetic surfaces. Since our calculation is not spin-polarized, the total charge density will always be shown.
import py4vasp
calc = py4vasp.Calculation.from_path( "./e10_STM_of_graphene/partial_charges/")
calc.partial_density.to_stm(selection="constant_current", current=8, supercell=[5,4])
We see the hexagonal symmetry of the graphene sheet and since both atoms in the system have the same chemical environment the tip height over both on-top positions is equivalent.
10.5 Questions¶
- What is the default mode for STM simulations in a py4vasp?
- If the current is increased from the default of 1nA, will the tip move closer to the surface or further away?
- Would it be possible to resolve spin-specific features of a surface, comparable to spin-polarized STM with py4vasp?