4.1. Nitrogen doped graphene (NDG)

Nitrogen-doped graphene (NDG) is a graphene sheet in which selected carbon sites are occupied or coordinated by N, tuning the local π-electron distribution and chemical reactivity. NDG combines high electrical conductivity, large specific surface area, and enhanced hydrophilicity, making it a useful platform for catalysis, sensing, and energy storage [Zhang17]. Depending on the different N configuration, distinct local density-of-states signatures that can be resolved by scanning tunnelling microscopy (STM) [Bar61, TH83, TH85].

In this tutorial we show how to use RESCU to simulate STM calculation based on Tersoff-Hamann methode. We will generate constant-height and constant-current images, extract one-dimensional STM line-scan profiles, and compute differential-conductance (\(dI/dV\)) maps, starting from the atomic structure shown in Fig. 4.1.1

../../../../_images/ndg.png

Fig. 4.1.1 Top view of N-doped graphene (graphitic substitution).

The calculation steps are given as follows:

  1. Prepare the N-doped graphene model (NDG);

  2. Self-consistent Field (SCF) calculation;

  3. Electronic analysis;

  4. STM simulations (Tersoff–Hamann);

    • Constant-height maps

    • Constant-current topography

  5. Line-scan profiles;

  6. Scanning Tunneling Spectroscopy (STS).

This NG model provides a compact workflow to obtain publication-quality STM constant-height and constant-current images, line-scan profiles, and \(dI/dV\) maps.

4.1.1. Build the structure

Generating the nitrogen doped graphene

The user can automatically generate the system by using the provided Python script gen_ndg.py. To execute this script, the user should have the ASE module properly installed.

Note

It is highly recommended to use an isolated environment, such as a Python virtual environment or a Conda environment, when installing additional Python modules on your system, in order to prevent potential conflicts with dependencies from other projects.

Generate the NDG structure by running:

python gen_ndg.py

A file should be generated (download here: ndg.cif). The structure is constructed using lattice parameters of a = b = 2.46 Å; c = 20 Å. The user can open this cif file using ase_gui or other viewer softare, like VESTA [KI2011].

4.1.2. Generate the SCF input file

To begin we’ll convert the ndg.cif file into a RESCU input file using the provided Python script cif2RESCU.py, as shown bellow:

python cif2RESCU.py -i ndg.cif

After running this command, two files should be generated:

These files will be used for the single-point calculation step. Go to the folder where the input files are located.

ndg.xyz file

50

C       0.000000000000000      0.000000000000000     10.000000000000000
C       1.230000000000000      0.710140831103240     10.000000000000000
C      -1.230000000000000      2.130422493309719     10.000000000000000
C      -0.000000000000001      2.840563324412958     10.000000000000000
C      -2.460000000000001      4.260844986619437     10.000000000000000
C      -1.230000000000001      4.970985817722677     10.000000000000000
C       8.609999999999999      6.391267479929156     10.000000000000000
C      -2.460000000000002      7.101408311032396     10.000000000000000
C      -4.920000000000002      8.521689973238875     10.000000000000000
C      -3.690000000000002      9.231830804342115     10.000000000000000
C       2.460000000000000      0.000000000000000     10.000000000000000
C       3.690000000000000      0.710140831103240     10.000000000000000
C       1.230000000000000      2.130422493309719     10.000000000000000
C       2.459999999999999      2.840563324412958     10.000000000000000
C      -0.000000000000001      4.260844986619437     10.000000000000000
C       1.229999999999998      4.970985817722677     10.000000000000000
C      -1.230000000000002      6.391267479929156     10.000000000000000
C      -0.000000000000002      7.101408311032396     10.000000000000000
C      -2.460000000000002      8.521689973238875     10.000000000000000
C      -1.230000000000002      9.231830804342115     10.000000000000000
C       4.920000000000000      0.000000000000000     10.000000000000000
C       6.150000000000000      0.710140831103240     10.000000000000000
C       3.690000000000000      2.130422493309719     10.000000000000000
C       4.920000000000000      2.840563324412958     10.000000000000000
N       2.459999999999999      4.260844986619437     10.000000000000000
C       3.690000000000000      4.970985817722677     10.000000000000000
C       1.229999999999998      6.391267479929156     10.000000000000000
C       2.459999999999999      7.101408311032396     10.000000000000000
C      -0.000000000000002      8.521689973238875     10.000000000000000
C       1.229999999999998      9.231830804342115     10.000000000000000
C       7.380000000000000      0.000000000000000     10.000000000000000
C       8.609999999999999      0.710140831103240     10.000000000000000
C       6.149999999999999      2.130422493309719     10.000000000000000
C       7.380000000000001      2.840563324412958     10.000000000000000
C       4.919999999999999      4.260844986619437     10.000000000000000
C       6.149999999999999      4.970985817722677     10.000000000000000
C       3.689999999999998      6.391267479929156     10.000000000000000
C       4.919999999999998      7.101408311032396     10.000000000000000
C       2.459999999999998      8.521689973238875     10.000000000000000
C       3.689999999999997      9.231830804342115     10.000000000000000
C       9.840000000000000      0.000000000000000     10.000000000000000
C      11.070000000000000      0.710140831103240     10.000000000000000
C       8.609999999999999      2.130422493309719     10.000000000000000
C       9.840000000000000      2.840563324412958     10.000000000000000
C       7.379999999999999      4.260844986619437     10.000000000000000
C       8.609999999999998      4.970985817722677     10.000000000000000
C       6.149999999999999      6.391267479929156     10.000000000000000
C       7.379999999999999      7.101408311032396     10.000000000000000
C       4.919999999999998      8.521689973238875     10.000000000000000
C       6.149999999999998      9.231830804342115     10.000000000000000

Note

The ndg.xyz file contains a total of 50 atoms, which the atom number 24 represent the N atom doping.

The scf input file

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[   12.300000     0.000000     0.000000        ];
[   -6.150000    10.652112     0.000000        ];
[    0.000000     0.000000    20.000000        ];
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Functional and LibXC Information
# -------------------------------------------------------------
functional.libxc = true % required for vdW
functional.list  = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}

# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = './pseudos/C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = './pseudos/N_ONCV_PBE_DZP.mat'
atom.xyz = 'ndg.xyz'

# ==============================================================
# ===================== SCF  ===================================
# ==============================================================

# -------------------------------------------------------------
# LCAO (Linear Combination of Atomic Orbitals) Configuration
# -------------------------------------------------------------
LCAO.status   = true
LCAO.mulliken = true

# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'self-consistent'
info.savepath        = 'results/ndg_scf'
info.outfile         = 'ndg_scf.out'

# -------------------------------------------------------------
# K-Point Grid and Smearing
# -------------------------------------------------------------
kpoint.type  = 'MonkhorstPack'
kpoint.gridn = [9,9,1]
kpoint.sigma = 0.05

# -------------------------------------------------------------
# Mixing Method for SCF Convergence
# -------------------------------------------------------------
mixing.method = 'broyden'
mixing.beta   = 0.05
mixing.type   = 'density'
mixing.tol    = [1e-5 1e-5]

# -------------------------------------------------------------
# Eigensolver Settings
# -------------------------------------------------------------
eigensolver.algo      = 'cfsi'
eigensolver.tol       = 5e-2
eigensolver.emptyBand = 36

# -------------------------------------------------------------
# Force Calculation Settings
# -------------------------------------------------------------
force.status = true

Now, we are ready to start the calculation. However, before doing that we’ll briefly present the simulation parameters specified in the generated file. The system information and simulation parameters are defined by keywords. These keywords generally have self-explanatory names, but you may refer to the input reference for more information.

The keywords specify the following:

  • option.timeLimit This keyword let RESCU knows how long (in hours) it should run before quitting

  • option.maxSCFiteration Maximal number of self-consistent steps; the Kohn-Sham solver will return an answer even if it has not converged.

  • option.precision This keyword sets various other keywords related to precision (e.g. domain.fourierInit, domain.highres, domain.vnlReal, diffop.method, interpolation.method, interpolation.vloc). It is not relevant for atomic orbital calculations (LCAO.status = true).

  • units.length Unit of length used for all structural inputs (atomic positions, lattice vectors, cell size). Accepted values are bohr (default) and angstrom.

  • units.energy Defines the default units of length in the input file.

  • domain.latvec Lattice vectors. The matrix must contain three row-vectors.

  • domain.boundary Determine the boundary conditions along each dimension.

  • domain.lowres controls the real space resolution.

  • functional.libxc Use the LibXC wrapper when functional.libxc = true to include functional from the LibXC.

  • functional.list List of exchange-correlation functional used in the calculation.

  • element(1).species It is a tag for the element structure. It is usually the chemical symbol of the corresponding element.

  • element(1).path is a string containing the path to the pseudo-potential file associated with the element.

  • atom.xyz Is an array of row-vectors containing the position of each atom.

  • LCAO.status If LCAO.status = true, RESCU performs a calculation using the atomic orbital basis. This is generally faster than performing a purely real space calculation, but the atomic orbital basis must be validated for correctness.

  • LCAO.mulliken Computes the Mulliken population matrix after converging the ground state density. An atomic orbital basis must be found in the pseudopotential file. If the wavefunctions are in real-space form, they will be projected onto the atomic orbital basis.

  • info.calculationType This input parameter specifies the type of the calculation that RESCU should perform.

  • info.savepath The path and name prefix of the file where the results will be writer. RESCU appends the rank number to the file name prefix, which for serial calculations is always 0.

  • info.outfile File in which writes the simulation information.

  • kpoint.type MonkhorstPack is used to perform integrals over the Brillouin zone; line is used to plot band structures.

  • kpoint.gridn sets the number of k-points used to sample the first Brillouin zone.

  • kpoint.sigma Use Fermi-Dirac smearing if option.sigma \(>\) 0; sigma in Hartree by default.

  • mixing.method If equal to linear and pulay the mixing method of the same name will be used. If equal to broyden method of Srivastava will be condidered.

  • mixing.beta Mixing fraction. This is used by all mixers.

  • mixing.type Determine whether to mix the density, density-matrix or the potential.

  • mixing.tol For self-consistent (ground state) calculations, the self-consistent procedure stops if the charge variation per valence electron is less than mixing.tol(1) and the total energy variation per valence electron is less than mixing.tol(2). Type [1,2].

  • eigensolver.algo Determine the algorithm used for real space diagonalization.

  • eigensolver.tol Maximal tolerance for the eigenvector residuals. When doing a self-consistent calculation, the eigensolvers use the minimum of eigensolver.tol and dRho/10.

  • eigensolver.emptyBand Number of bands included in the calculation in addition to the number of valence bands

  • force.status Determines whether to compute the forces. It is sometimes costly to compute forces, and hence one may choose to avoid it. The default is true for real space self-consistent calculations and false otherwise.

Now, we are ready to start the calculations. Copy the C and N Pseudopotential file from the database to ./pseudos directory. Execute RESCU and use the -i option to point to the input file.

rescu -i ndg_scf.input

After the calculation is finished, the following output files are generated:

  • ndg_scf.out;

  • .results/ndg_scf.h5;

  • .results/ndg_scf.mat.

The file ndg_scf.mat contains the converged ground state density that will be necessary to compute the PDOS.

4.1.3. Electronic analysis

Now we’ll investigate the electronic states from ndg by calculatinting the pdos with the input below.

The pdos input file

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[      6.153000    -10.657309      0.000000        ];
[      6.153000     10.657309      0.000000        ];
[      0.000000      0.000000     20.000000        ]
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Funtional and LiXC Information
# -------------------------------------------------------------

#functional.libxc      = true % required for vdw
functional.list      = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}


# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = 'C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = 'N_ONCV_PBE_DZP.mat'
atom.xyz             =  'grafN.xyz'

# ==============================================================
# ===================== PDOS ===================================
# ==============================================================

# -------------------------------------------------------------
# PDOS Calculation Input
# -------------------------------------------------------------
rho.in           = 'results/ndg_scf'
dos.status       = true
dos.range        = [-1, 1]
dos.resolution   = 0.01
dos.projAtom     = 24


# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'dos'
info.savepath        = 'results/ndg_pdos'
info.outfile         = 'ndg_pdos.out'


# -------------------------------------------------------------
# K-Point Grid and Smearing
# -------------------------------------------------------------
kpoint.type  = 'MonkhorstPack'
kpoint.gridn = [9,9,1]
kpoint.sigma = 0.005

The PDOS input file intentionally repeats the structural section (units, lattice vectors, boundary conditions, elements, and atom.xyz). Although rho.in points to the converged density from the SCF step, RESCU still rebuilds the Hamiltonian and the PDOS projectors from the geometry and atomic basis defined in the current input. Re-stating these parameters makes the PDOS run self-contained and ensures it uses exactly the same system as the SCF calculation. Compared with SCF, we only adjust the analysis settings: switch info.calculationType = dos, dos.status = true, choose an energy window and resolution, and often tighten the smearing (e.g., kpoint.sigma = 0.005) to sharpen spectral features while keeping the same k-point mesh.

The specific keywords related to the PDOS analysis are given by:

  • dos.status Determines whether to compute the density of states.

  • dos.range 1 x 2 array giving the energy range, relative to the Fermi energy, in Hartree, for which the DOS is calculated.

  • dos.resolution Determines the resolution of the DOS energy axis; in Hartree by default.

  • dos.projAtom Calculate the density of states projected on the atomic orbitals of specific atoms.

Execute RESCU by:

rescu -i ndg_pdos.input

After the calculation finish the user can quicly plot the results with:

rescu -p .results/ndg_pdos.mat

PDOS analysis

The PDOS .fig file generated is given by:

../../../../_images/ndg_pdos_plot_pdos.png

Fig. 4.1.2 Projected Density of states of GDN.

As one can see, the density of states is overwhelmingly dominated by C-projected states (black), while the N contribution (red) remains small across the entire window. The finite DOS exactly at \(E_{F}\) and the strong C-derived resonance just below indicate metallic or degenerately doped behavior [Chen14].

The PDOS supports a picture of N-doped graphene where the dopant perturbs but does not localize states near \(E_{F}\), producing an n-type shift and preserving the delocalized \(\pi\) character of the electronic structure. No sharp, isolated N-localized impurity level appears in the gap region; instead, N states weakly hybridize with the carbon \(\pi\) network.

The ldos input file

For the subsequent STM analysis, we have now to re-do the PDOS calculation to generate the LDOS data. To do this make a copy of ndg_pdos.input to ndg_ldos.input. Comment the line given by dos.projAtom = 24 and add the following block:

## -------------------------------------------------------------
## LDOS Calculation Settings
## -------------------------------------------------------------
dos.ldos = true
dos.ldosShiftEf = true
##dos.ldosERange = [0,0.5]
dos.ldosEnergy  = [-2.0:0.05:2.0]

The LDOS input file can download here: ndg_ldos.input.

The specific keywords related to the LDOS analysis are given by:

  • dos.ldos Determines whether to compute the density of states.

  • dos.ldosShiftEf If dos.ldosShiftEf = true, then the values given in dos.ldosEnergy are relative to the Fermi energy.

  • dos.ldosERange The ldos is integrated over the energy range given by dos.ldosERange. This is also refered to as the partial charge. The range is relative to the energy given by the keyword dos.ldosEnergy, whose default value is the Fermi energy. Note that dos.ldosEnergy is itself relative to the Fermi energy if dos.ldosShiftEf = true.

  • dos.ldosEnergy Energies at which the local density of states is calculated. If dos.ldosEnergy (or alternate keywords) is not specified, the LDOS is calculated at the Fermi energy. Note that this will typically lead to a zero LDOS in semiconductors. The energies are relative to the Fermi energy if dos.ldosShiftEf = true. The Fermi energy is recalculated during the DOS calculation, and hence it may differ from the result of the self-consistent calculation. This keyword has precedence over dos.ldosEnergyi and dos.ldosEnergyk.

Note

This PDOS calculation also generated LDOS on discrete energy channels from −2 to +2 eV (dos.ldosEnergy = [-2:0.05:2]), which is the range where dopant-related features are observed.

We will use dos.ldosEnergy to specify discrete energy channels for the later \(dI/dV\) step. This produces energy-resolved LDOS data that can be memory-intensive, since storage scales with the number of energies, k-points, bands, and real-space grid points. As a lighter alternative, you may use dos.ldosERange, which integrates all states within a single energy window into one channel. All subsequent STM/LDOS analyses that rely on integrated density will work the same way, but \(dI/dV\) cannot be plotted from a single integrated channel because the energy resolution is lost.

Execute RESCU by:

rescu -i ndg_ldos.input

4.1.4. STM simulations maps

The LDOS computed in the previous step is saved to ./results/ndg_ldos.mat and ./results/ndg_ldos.h5. These files contain the spatially resolved LDOS on the real-space grid and are used to generate constant-height and constant-current STM maps within the Tersoff-Hamann approximation [TH83, TH85].

4.1.4.1. Constant-height maps

First the we’ll get the averaged current at z = 14.4 Å. (the graphene layer is at z = 10.0 Å). To compute constant-height STM maps, use the input file below.

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[   12.300000     0.000000     0.000000        ];
[   -6.150000    10.652112     0.000000        ];
[    0.000000     0.000000    20.000000        ];
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Functional and LibXC Information
# -------------------------------------------------------------
functional.libxc = true % required for vdW
functional.list  = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}

# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = './pseudos/C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = './pseudos/N_ONCV_PBE_DZP.mat'
atom.xyz = 'ndg.xyz'

# ==============================================================
# ===================== STM-TERSOFF ============================
# ==============================================================

# -------------------------------------------------------------
# STM Calculation Input
# -------------------------------------------------------------
stm.ldos.in         = 'results/ndg_ldos'

# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'stm-tersoff'
info.savepath        = 'results/ndg_stmh'
info.outfile         = 'ndg_stmh.out'

# -------------------------------------------------------------
# Scan mode
# -------------------------------------------------------------
stm.scanMode    = 'H'

# -------------------------------------------------------------
# Height mode selection
# -------------------------------------------------------------
stm.heightMode  = 1
stm.h           = 14.42

# -------------------------------------------------------------
# Geometric offset
# -------------------------------------------------------------
stm.zOffset     = 0

# -------------------------------------------------------------
# Contrast / percentile
# -------------------------------------------------------------
stm.p           = 20

# -------------------------------------------------------------
# Bias integration
# -------------------------------------------------------------
stm.bias.mode   = 'window'
stm.bias.V      = -0.4

# -------------------------------------------------------------
# Tiling and aesthetics
# -------------------------------------------------------------
stm.m           = 3
stm.n           = 3
stm.colormapId  = 5
stm.outPrefix   = 'stm-N04'

Tip

As with the PDOS input, we embed the NDG structural data directly in the file so the example remains fully self-contained and reproducible, preserving all geometry and calculation parameters without relying on external files.

The STM input file for constant-height can download here: ndg_stm-constH.input.

The specific keywords related to the STM constant-height analysis are given by:

  • stm.scanMode Determines the type of scan mode. H- constant‑height: fixed tip height; returns a current map (or raw LDOS). C - constant‑current: adjusts height to reach a current/LDOS setpoint. S - 2D slice: renders a 2D LDOS slice at a fixed height.

  • stm.heightMode Specifies how the tip height is provided: 1 - single height stm.h. 2 - uniform scan from stm.h1 to stm.h2 with stm.npoints (inclusive).

  • stm.h Tip height in Å above the reference plane stm.zOffset when stm.heightMode = 1.

  • stm.zOffset Geometric reference plane in Å (e.g., topmost atomic layer).

  • stm.p Percentile used to define the setpoint when stm.setpoint = 'percentile' (in constant‑current mode).

  • stm.bias.mode Controls the energy integration of the LDOS: window or none.

  • stm.bias.V Sample bias in eV (negative for occupied‑state imaging).

  • stm.m and stm.n Periodic tiling factors along lattice vectors a (m) and b (n) applied before rendering STM images.

  • stm.colormapId Colormap selector for rendering

  • stm.outPrefix Base prefix for output file names (images and dumps).

Now execute RESCU by:

rescu -i ndg_stm-constH.input

After the calculation finishes, the following output files are generated:

  • ndg_stmh.out;

  • stm-N04_H_14.420.png;

STM map analysis

After running the calculation, the STM constant-height map shows the simulated tunneling current over the NDG surface for a tip held at a fixed height of 14.4 Å (4.4 Å above the sheet). Note that we use this value of heigth in our simulation to get a map that is in line with experimental result [Zhao11]. The x and y axes are the in-plane atomic coordinates (in Å), and the color scale on the right encodes the magnitude of the current I in arbitrary units

../../../../_images/stm-N04_H_14.420.png

Fig. 4.1.3 The STM constant-height map.

Brighter (yellow/white) spots correspond to regions of higher local density of occupied states within the chosen bias window [-0.4, 0] eV relative to the Fermi level, meaning these sites contribute more strongly to the tunneling current at this bias [Joucken12]. These high-intensity spots are associated with the local electronic states around the dopant environment, while the darker regions correspond to areas with lower LDOS in that energy range.

Tip

As an additional exercise, users can also generate the map at a bias of +0.4 V. To do so, simply set:

stm.bias.V      =  0.4
stm.outPrefix   = 'stm-P04'

4.1.4.2. Constant-current maps

To compute constant-current STM maps the users can follow the previous step using the input file below. In this case, to obtain a simulated constant-current topography comparable to an experimental STM image [Joucken12], we set stm.h = 12.80 (in Å), which defines the reference tip height for the feedback loop.

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[   12.300000     0.000000     0.000000        ];
[   -6.150000    10.652112     0.000000        ];
[    0.000000     0.000000    20.000000        ];
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Functional and LibXC Information
# -------------------------------------------------------------
functional.libxc = true % required for vdW
functional.list  = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}

# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = './pseudos/C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = './pseudos/N_ONCV_PBE_DZP.mat'
atom.xyz = 'ndg.xyz'


# ==============================================================
# ===================== STM-TERSOFF ============================
# ==============================================================

# -------------------------------------------------------------
# STM Calculation Input
# -------------------------------------------------------------
stm.ldos.in         = 'results/ndg_ldos'

# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'stm-tersoff'
info.savepath        = 'results/ndg_stm'
info.outfile         = 'ndg_stm.out'

# -------------------------------------------------------------
# Scan mode
# -------------------------------------------------------------
stm.scanMode    = 'C'

# -------------------------------------------------------------
# Height mode selection
# -------------------------------------------------------------
stm.heightMode   = 1
stm.h            = 12.8

# -------------------------------------------------------------
# Geometric offset
# -------------------------------------------------------------
stm.zOffset     = 0

# -------------------------------------------------------------
# Contrast / percentile (used in cnt-I setpoint='percentile')
# -------------------------------------------------------------
stm.p           = 20

# -------------------------------------------------------------
# Bias integration (requires LDOS with energy axis)
# -------------------------------------------------------------
stm.bias.mode   = 'window'
stm.bias.V      = -0.4

# -------------------------------------------------------------
# Tiling and aesthetics
# -------------------------------------------------------------
# stm.images.enabled = 1
stm.m           = 3
stm.n           = 3
stm.colormapId  = 3
stm.outPrefix   = 'stmC-N04'

Now execute RESCU by:

rescu -i ndg_stm-constC.input

After the calculation finishes, the following output files are generated:

  • ndg_stmc.out;

  • stmC-N04_C_12.800.png;

Note

The user may scan a range of tip heights along the z direction to identify a suitable operating distance for the STM tip. This is done by enabling height sweep mode and including the following keywords:

'stm.heightMode = 2'

'stm.h1'   (Lower  tip heights in Å)

'stm.h2'   (Upper tip heights in Å)

'stm.npoints'   (Number of heights (including endpoints) in the range)

STM map analysis

After running the STM constant-current calculation, the user should obtain a topographic map shown below. In this image, the x and y axes (in Å) correspond to the in-plane atomic coordinates of the simulated surface.

../../../../_images/stmC-N04_C_12.800.png

Fig. 4.1.4 The STM constant-height map.

The color scale represents the apparent tip height z (in Å) required to maintain a fixed tunneling current (“setpoint=auto”) within the chosen bias window [-0.4,0] eV. Warmer colors (yellow/red) indicate regions where the tip must retract to keep the current constant, these appear as bright protrusions in STM. Cooler colors (blue) indicate lower apparent height, where the tip can approach the surface more closely.

The user should see pronounced triangular features at specific sites. These correspond to electronically enhanced regions (e.g. around the dopant environment), which appear higher in constant-current mode. The surrounding honeycomb pattern reflects the graphene-like lattice contrast of the host layer.

4.1.5. STM linescans

An important topography analises from STM data is given by a line scan, a 1D cut through the STM data. Instead of looking at the full 2D map, we’ll select a path (a line) across the surface and record how the signal changes along that path.

Tip

Why this is useful:

  • It gives quantitative contrast between two features that look “bright” or “dim” in the map;

  • It lets you measure lateral distances (e.g. defect–defect spacing, feature size).

  • It helps identify whether an observed protrusion is really topographic (geometric height) or electronic (LDOS enhancement).

To compute linescans the users can follow input file below.

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[   12.300000     0.000000     0.000000        ];
[   -6.150000    10.652112     0.000000        ];
[    0.000000     0.000000    20.000000        ];
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Functional and LibXC Information
# -------------------------------------------------------------
functional.libxc = true % required for vdW
functional.list  = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}

# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = './pseudos/C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = './pseudos/N_ONCV_PBE_DZP.mat'
atom.xyz             = 'ndg.xyz'

# ==============================================================
# ===================== STM-TERSOFF ============================
# ==============================================================

# -------------------------------------------------------------
# STM Calculation Input
# -------------------------------------------------------------
stm.ldos.in         = 'results/ndg_ldos'

# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'stm-tersoff'
info.savepath        = 'results/ndg_stm-ls'
info.outfile         = 'ndg_stmls.out'

# -------------------------------------------------------------
# Scan mode
# -------------------------------------------------------------
stm.scanMode               = 'C'

# -------------------------------------------------------------
# Height mode selection
# -------------------------------------------------------------
stm.heightMode  = 1
stm.h           = 12.80

# -------------------------------------------------------------
# Geometric offset
# -------------------------------------------------------------
stm.zOffset     = 0

# -------------------------------------------------------------
# Contrast / percentile (used in cnt-I setpoint='percentile')
# -------------------------------------------------------------
stm.p           = 20

# -------------------------------------------------------------
# Bias integration (requires LDOS with energy axis)
# -------------------------------------------------------------
stm.bias.mode   = 'window'
stm.bias.V      = -0.4

# -------------------------------------------------------------
# Tiling and aesthetics
# -------------------------------------------------------------
stm.m           = 3
stm.n           = 3
stm.colormapId  = 3
stm.outPrefix   = 'stmls'

# -------------------------------------------------------------
# Line-scan (extracted along a segment on the rendered map)
# -------------------------------------------------------------
stm.linescan.enabled = 1
stm.linescan.n      = 600
stm.linescan.x1     = 0
stm.linescan.y1     = 0
stm.linescan.x2     = 18.45000
stm.linescan.y2     = 31.95634
stm.linescan.width  = 0.0

The line scan input file can be download here: ndg_stm-ls.input. The specific keywords related to the linescans analysis are given by:

  • stm.linescan.enabled Enables line-scan extraction after rendering the STM map.

  • stm.linescan.n Number of uniformly spaced samples taken along the line segment from (x1, y1) to (x2, y2) in the plotting plane.

  • stm.linescan.(x1, x2), stm.linescan.(y1, y2) Start point (x1, y1) and end point (x2, y2) of the line in Å, expressed in the STM image plane (after tiling by stm.m, stm.n and at the effective height given by stm.zOffset + h).

  • stm.linescan.width Transverse averaging width in Å around the line (perpendicular direction). If 0, no lateral averaging is performed. Small values help reduce noise by averaging pixels within a narrow strip.

To calculate the linescans execute RESCU by:

rescu -i ndg_stm-ls.input

After the calculation finishes, the following output files are generated:

  • ndg_stm-ls.out;

  • stmls_C_12.800.png;

  • stmls_C_12.800_linescan.png.

Linescans analysis

This is the constant-current STM topography map (same as before), but now it includes a dashed line connecting two points (user-defined): the starting point (marked with a circle) and the ending point (marked with a square). This dashed line is the path along which the line scan is extracted.

../../../../_images/stmls_C_12.800.png

Fig. 4.1.5 Dashed line is the path along which the line scan is extracted.

This plot shows the apparent tip height z (in Å) as a function of the distance s along the dashed path (in Å). Each tall peak corresponds to the tip moving upward over a bright protrusion in the STM image [Zhao11]. The flatter baseline between peaks reflects the lower background surface.

The periodic spacing between the peaks gives you a direct measurement of the separation between equivalent features (for example, between dopant sites) in real space.

../../../../_images/stmls_C_12.800_linescan.png

Fig. 4.1.6 The STM line-scan.

In practice, this line-scan output can be compare directly to an experimental STM height trace along a chosen cut on the surface.

4.1.6. Scanning tunneling spectroscopy

Scanning tunneling spectroscopy (STS) is an extension of scanning tunneling microscopy (STM) that measures not only how the surface “looks” in real space, but also how its electronic states are distributed in energy. In practice, STS records the tunneling current (\(I\)) as a function of the applied sample bias (V) at a fixed tip position. By taking the numerical derivative \(dI/dV\), we obtain a signal that is approximately proportional to the local density of states (LDOS) of the sample at that position and at the corresponding energy \(E = E_{F} + eV\)

In this workflow, we evaluate the LDOS on a grid above the surface and use it to predict the differential conductance signal. For this NDG system, we’ll evaluate the STS signal between -1.5 eV to +1.5 eV (Fig. 4.1.2)

To compute the STS the users can follow input file below.

# -------------------------------------------------------------
# Calculation Options
# -------------------------------------------------------------
option.timeLimit       = 2
option.maxSCFiteration = 400
option.precision       = 'real-med'

# -------------------------------------------------------------
# Unit Settings
# -------------------------------------------------------------
units.length = 'Angstrom'
units.energy = 'eV'

# -------------------------------------------------------------
# Domain Settings (Lattice Vectors and Boundaries)
# -------------------------------------------------------------
domain.latvec = [
[   12.300000     0.000000     0.000000        ];
[   -6.150000    10.652112     0.000000        ];
[    0.000000     0.000000    20.000000        ];
]
domain.boundary = [1,1,2]
domain.lowres   = 0.3

# -------------------------------------------------------------
# Functional and LibXC Information
# -------------------------------------------------------------
functional.libxc = true % required for vdW
functional.list  = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}

# -------------------------------------------------------------
# Atomic and Elemental Information
# -------------------------------------------------------------
element(1).species   = 'C'
element(1).path      = './pseudos/C_ONCV_PBE_DZP.mat'
element(2).species   = 'N'
element(2).path      = './pseudos/N_ONCV_PBE_DZP.mat'
atom.xyz = 'ndg.xyz'

# ==============================================================
# ===================== STM-TERSOFF ============================
# ==============================================================

# -------------------------------------------------------------
# STM Calculation Input
# -------------------------------------------------------------
stm.ldos.in         = 'results/ndg_ldos'

# -------------------------------------------------------------
# General Calculation Information
# -------------------------------------------------------------
info.calculationType = 'stm-tersoff'
info.savepath        = 'results/ndg_sts'
info.outfile         = 'ndg_sts.out'

# -------------------------------------------------------------
# Bias integration (requires LDOS with energy axis)
# -------------------------------------------------------------
stm.bias.mode   = 'dIdV'
stm.sts.enabled   = 1
stm.sts.x         = 0
stm.sts.y         = 0
stm.sts.z         = 14.44
stm.sts.bias0     =  -1.5
stm.sts.bias1     =  1.5
stm.sts.biasstep  =  0.05

The STS input file can download here: ndg_STS.input.

The specific keywords related to the STS analysis are given by:

  • stm.bias.mode Controls the energy integration of the LDOS emulating a sample bias: ‘window’ — integrates from the Fermi level to V (or from V to the Fermi level if V < 0). ‘none’ — disables energy integration; select a fixed LDOS channel with stm.ldosChannel.

  • stm.sts.enabled Enables the generation of image outputs (PNG/MAT, etc.). Set to 0 to run numerical/STS only without rendering figures.

  • stm.sts.x, stm.sts.y Absolute Å coordinates on the rendered STM image plane at which the STS spectrum is sampled (after tiling by stm.m, stm.n). Not fractional coordinates.

  • stm.sts.z Absolute height in Å at which the spectrum is taken. If omitted, defaults to the effective tip height zOffset + h (when available from the current scan settings).

  • stm.sts.bias0, stm.sts.bias1, and stm.sts.biasstep Bias sweep definition in eV. The spectrum is sampled from bias0 to bias1 (inclusive when possible) using uniform increments biasstep. Negative bias probes occupied states; positive bias probes empty states.

Now the user can execute RESCU by:

rescu -i ndg_STS.input
After running the STS calculation, the user will obtain the plots:

  • ndg_sts.out;

  • stm_STS_x0.000_y0.000_z14.440_IV.png;

  • stm_STS_x0.000_y0.000_z14.440_dIdV.png.

STS analysis

The Fig. 4.1.7 shows the tunneling current \(I\) (in nA) as a function of the bias \(V\) (in eV), at the same spatial position. At negative and small positive bias, the current remains close to zero, meaning that the tip-sample junction is effectively “off” because there are not many accessible states in that energy range. Beyond ~0.5 to 1.0 eV, the current starts to rise and then increases rapidly.

../../../../_images/stmSTS-IV.png

Fig. 4.1.7 STS: \(I-V\) spectrum.

The Fig. 4.1.8 shows the differential conductance \(dI/dV\) (in nA/eV) as a function of the applied sample bias \(V\) (in eV),, evaluated at a single tip position above the surface (x,y,z) = (0.00,0.00,14.44 Å). In the Tersoff-Hamann picture, \(dI/dV\) is proportional to the local density of states (LDOS) of the sample at that point in space and energy. The regions where \(dI/dV ≈ 0\), indicating very few available states to tunnel into/out of, and sharp onsets or peaks at positive bias.

The rapid increase in \(dI/dV\) above ~1 eV reflects the emergence of unoccupied states that strongly contribute to tunneling. These features can be directly compared to experimental STS curves measured over a specific atomic site or defect.

../../../../_images/stm_STS_dIdV.png

Fig. 4.1.8 STS: \(dI/dV\) spectrum.

Together, these two plots give both the integral view (\(I-V\)) and the energy-resolved view (\(dI/dV\)) of the local electronic structure at the chosen point above the surface. This is the simulated STS signal you would compare directly to experimental spectroscopy taken with the STM tip parked over the same type of site.

4.1.7. References

[Zhang17]

Y. Q. Han, Y. L. Zhang, X. Zhou, B. Deng, Q. Li, M. Liu, J. Zhao, Z. Liu, and O. Zhang, Strong Adlayer–Substrate Interactions “Break” the Patching Growth of h-BN onto Graphene on Re(0001) ACS Nano 11, (2017), p. 1807.

[Bar61]

J. Bardeen, Tunneling from a Many-Particle Point of View Phys. Rev. Lett. 6, (1961), p. 57.

[TH83]

J. Tersoff and D. R. Hamann, Theory of the scanning tunneling microscope Phys. Rev. Lett. 50, (1983), p. 1998.

[TH85]

J. Tersoff and D. R. Hamann, Theory of the scanning tunneling microscope Phys. Rev. B 31, (1985), p. 805.

[KI2011]

K. Momma and I. Fujio, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data Journal of Applied Crystallography 44, (2011), p. 1272.

[Chen14]

X. Wang, G. Sun, P. Routh, D.-H .Kim, W. Huangb and P. Chen, Heteroatom-doped graphene materials: syntheses, properties and applications Chem Soc Rev 43, (2014), p. 7067.

[Joucken12] (1,2)

F. Joucken et al., Localized state and charge transfer in nitrogen-doped graphene Phys. Rev. B 85, (2012), p. 161408.

[Zhao11] (1,2)

L. Zhao et al., Visualizing Individual Nitrogen Dopants in Monolayer Graphene Science 333, (2011), p. 999.