Skip to main content

Bulk and Interface Transport

The NEGF + TBtrans pipeline is not reserved for molecular junctions. As long as the electrode/scattering-region layout is maintained, perfect crystals, metal interfaces, and 2D sheets can all be computed with the same workflow. This chapter organizes three types of general, molecule-free systems, points out that — in exact opposition to the 1D vacuum system — the transverse k-average becomes part of the result, and maps out how far each chapter of the existing 1D tutorial carries over unchanged.

Learning objectives

  • Explain that in a perfect-crystal calculation, where electrode and device are the same material, T(E)T(E) becomes a counting of bands (channels).
  • Describe the device-construction principles of interface/heterojunction and 2D-sheet transport in the same language as the 1D junction.
  • Write the transverse k-average definition of T(E)T(E) and explain why 3D/2D systems are less sensitive to η\eta.
  • Apply the procedural checklist for starting transport with a new material or new electrode.

Three representative types

Type 1 — perfect crystal: the reference system for ballistic verification

If the electrode and the device are the same crystal of the same material, there is no scatterer. Then the transmission at each k\mathbf{k}_\parallel is exactly equal to the number of Bloch bands open along the transport direction at that energy.

T(E,k)=Nband(E,k){0,1,2,}T(E, \mathbf{k}_\parallel) = N_{\mathrm{band}}(E, \mathbf{k}_\parallel) \in \{0, 1, 2, \dots\}

That is, the T(E)T(E) of a perfect crystal is another way of counting the band structure (Chapter 02). You can sketch it in advance on the band diagram — at a fixed k\mathbf{k}_\parallel, draw a horizontal line at energy EE and count how many times it crosses the transport-direction bands (the number of states with positive propagation direction); that count is T(E,k)T(E, \mathbf{k}_\parallel).

(at fixed k∥) number of bands crossed by the energy horizontal line -> T(E, k∥)
region where 2 bands pass 2
region where only 1 band passes 1
gap (no bands) 0

Its significance as a calculation lies in setup verification. When starting transport with a new electrode or new material, the first calculation must always be the pristine system, and if T(E,k)T(E, \mathbf{k}_\parallel) fails to recover the integer staircase above, the problem is the setup (excessive η\eta, electrode–device lattice mismatch, principal layer violation), not physics. It is the generalization of the verification done with the integer plateaus of the pristine chain in the 1D example (Chapter 07). In bulk, the direction is a variable too — even for the same crystal, the number of channels open along the transport axis changes with the crystallographic direction chosen, so T(E)T(E) is a direction-dependent quantity.

Type 2 — interfaces and heterojunctions

A metal A/metal B junction, or a defect layer inside the same material (a stacking fault, an impurity layer), belongs to this type. Assemble the device as [A electrode replica] + [interface/defect region + buffer layers] + [B electrode replica], and have the left and right TS.Elec blocks point to different TSHS files — the block structure itself is the same as in Chapter 06. Two points to watch:

  • Transverse lattice matching. The surface cells of the two materials must share the same transverse lattice. If the lattices differ, one side has to be strained to match, and since the strain itself changes the electronic structure, you must state which side was strained and by how much.
  • A defect layer inside the same material is simpler. When the electrodes are identical on both sides, as for stacking faults and impurity layers, there is no lattice-matching problem and the layout is exactly that of the 1D example's C19N (pristine electrodes + central substitution) — just transfer the device-assembly logic of Chapter 06 into 3D.
  • Interface dipole and buffer length. When two metals with different work functions meet, charge transfer and a dipole appear at the interface. That electrostatic perturbation must finish decaying inside the device, so the bulk-matching judgment of the boundary-layer potential/DOS (scattering region convergence) applies unchanged.
The SCF resolves the alignment of the two electrodes' Fermi levels

At 0 V the two electrodes share a common chemical potential (μL=μR=EF\mu_L = \mu_R = E_F). Even if the work functions of materials A and B differ, the user need not align the levels by hand — the TranSIESTA SCF redistributes charge at the interface to create a dipole, and that dipole absorbs the electrostatic potential offset between the two bulks. What the user must check is not the alignment itself but whether the perturbation of that dipole decayed before reaching the device boundary.

Type 3 — surfaces and 2D sheets

In-plane transport in a 2D material such as graphene differs only in the arrangement of directions: the transport direction (periodic, set k = 1 and let the self-energy handle it), the in-plane transverse direction (periodic, k dense), and the out-of-plane direction (vacuum, k = 1). It is thus an intermediate configuration between the 1D chain (all transverse directions vacuum) and 3D bulk (all transverse directions periodic), and the transverse k-average is carried out over one direction only. In graphene, the crystallographic direction of the transport axis (zigzag/armchair) is also a design variable that must be stated, and since the states near EFE_F are concentrated around the Dirac point (a specific k\mathbf{k} of the BZ), a coarse transverse k misses the low-energy transmission entirely — the reason k convergence tests are especially important in 2D systems.

Summarizing the direction arrangements of the three types together with the 1D example:

Systemtransporttransverse periodicvacuumTBtrans k-grid form
1D chain (Chapters 05–09)zznonexx, yy1×1×11\times1\times1
2D sheet (type 3)xxyyzz1×n×11\times n\times1
3D bulk (types 1, 2)zzxx, yynonen×n×1n\times n\times1

Normalization — T is not directly comparable when cell sizes differ

Transmission is the number of channels passing through the transverse cell, so it grows in proportion to the cell size. It must be normalized per unit area in 3D bulk (T/AT/A) and per unit width in a 2D sheet (T/WT/W) before it can be compared with calculations or literature using a different transverse cell. If the calculation examines the scattering by a single defect, the ratio to pristine (T/TpristineT/T_{\mathrm{pristine}}) is conversely a metric less sensitive to cell size — state in the axis label which normalization was used.

The transverse k-average is part of the result

In a 1D vacuum system the k-points were merely a convergence parameter (a single Γ\Gamma sufficed), but in a system with transverse periodicity the observable is a k-average.

T(E)=1NkkT(E,k)T(E) = \frac{1}{N_{k}} \sum_{\mathbf{k}_\parallel} T(E, \mathbf{k}_\parallel)

Since the band structure differs for each k\mathbf{k}_\parallel, the integer staircase of the open-channel count T(E,k)T(E, \mathbf{k}_\parallel) steps up and down at different energies for each k, and averaging yields a smooth curve. Thanks to this natural averaging, 3D/2D systems are relatively less sensitive to the numerical broadening η\eta — the sharp structures of individual k (band edges, van Hove singularities) get washed into the average. Conversely, the 1D vacuum system had no k-average and was therefore especially sensitive to η\eta convergence (this is the same content as "The special nature of 1D systems" in the η broadening artifact page, seen from the opposite side). Being less sensitive is not exemption, so the η\eta series and k series convergence checks are still performed.

The quantity that decides k-convergence

Transverse k convergence is judged not by "has the total energy converged?" but by "does the T(E)T(E) curve overlap with the one from a denser grid?" Check with an overlay over the whole energy window of interest, not with the single point T(EF)T(E_F) — the same principle as the length-convergence judgment in scattering region convergence.

k-resolved transmission is extracted from trans.TBT.nc by exactly the sisl route used in Chapter 07. The k-average is the default, and the kavg argument returns the per-k curves.

Extracting k-resolved T(E) (sisl)
import numpy as np
from sisl.io.tbtrans import tbtncSileTBtrans

t = tbtncSileTBtrans("trans.TBT.nc")

E = t.E # E - E_F (eV)
T_avg = t.transmission(*t.elecs[:2]) # k-average (default)
T_k = t.transmission(*t.elecs[:2], kavg=False) # per k -- the integer staircase is visible
w = t.wkpt # k-point weights

# Self-check: the weighted average must match the k-average
assert np.allclose((w[:, None] * T_k).sum(axis=0), T_avg)

For a perfect crystal, each row of T_k must be an integer staircase and T_avg a smooth curve — checking both properties at once is the setup verification of type 1 above.

The cheapest system for practice — the Au monatomic chain

Before building heavy bulk example inputs, the feel of metallic electrodes can be acquired most cheaply with a 1D Au monatomic chain. The structure is the carbon chain of Chapter 05 with the element changed to Au (transverse vacuum, periodic along the transport direction), and the lattice constant is fixed by relaxing the chain itself following the procedure of Chapter 03. Since Au is an element for which relativistic effects enter the lattice and bands substantially, check that a scalar-relativistic pseudopotential is used. This system has only one band, 6s, crossing EFE_F, making it a single-channel system — a T(EF)=1T(E_F) = 1 plateau is expected in the pristine case, corresponding to the measured conductance of 1G01\,G_0 for real Au atomic wires. It is a good exercise in contrast with the doubly degenerate π\pi channels of the carbon chain (T(EF)=2T(E_F) = 2). Going below EFE_F, bands of 5d origin add channels. Since this system is still a transverse-vacuum 1D system, the Chapter 05–07 workflow runs unchanged without a k-average, making it a suitable stepping stone before moving to 3D Au electrodes (metal-electrode molecular junction).

Procedural checklist — starting transport with a new material

  1. Choose the electrode cell — fix the lattice by bulk relaxation and take the minimum periodic unit along the transport direction (for fcc [111], the ABC 3-layer period; an integer multiple if needed).
  2. Verify the principal layer — judge with the orbital reach whether hopping beyond adjacent cells is absent (Chapter 05 condition 2). If it is ambiguous, extend the cell along the transport direction.
  3. Assemble the device — electrode replica regions + region of interest (interface/defect) + buffer layers. Make the electrode replica regions atomically identical to the electrode cell, and set the buffer length by the criteria of scattering region convergence.
  4. Pristine verification — first confirm the integer structure of T(E)T(E) with a perfect crystal free of defects and interfaces (type 1). If it fails here, no subsequent result can be trusted.
  5. Transverse k convergence test — raise the TBtrans k stepwise (4×48×816×164\times4 \to 8\times8 \to 16\times16 and so on) and judge with a T(E)T(E) overlay. Place the TBT.kgrid_Monkhorst_Pack block at the very end of input.fdf (Chapter 07 pitfall 2).

Practical note — where the cost grows

Going from a 1D chain to 3D bulk changes the structure of the computational cost. Three factors compound multiplicatively.

  • Matrix size — the device Green function is cubically expensive in the number of scattering-region orbitals. In 3D there are many atoms per layer (9 per layer for 3×3 Au(111)), so the scattering region has far more orbitals than in 1D, which makes the saving of narrowing the Green-function solution region to the scattering region with TBT.Atoms.Device (Chapter 07) more important than in 1D.
  • Number of k-points — each transverse k is an independent calculation, so the cost is linear in NkN_k. On the other hand, there is no communication between k-points, so parallelization is easy (Chapter 12).
  • Number of energy points — linear in the energy grid, as in 1D. Estimate the total cost roughly as the product (energy points) × (k-points) × (matrix operations).

Estimating this product in advance when planning convergence tests (buffer length, k, η\eta) lets you judge, before computing, which axis's convergence series is expensive — and therefore which axis to fix first.

Mapping to the 1D tutorial chapters

Here is how far the procedures of the existing chapters remain valid in transverse-periodic systems.

ChapterApplicabilityRemarks
05 Electrodek-sampling only differsA dense transverse k is merely added; the principal layer judgment and the TS.HS.Save procedure are identical
06 Device 0 Vk-sampling only differsTransverse k in the device SCF (n×n×1n\times n\times1). The TS.ChemPots/TS.Elecs/contour block structure is identical
07 TBtransCautionThe minimal option set is the same, but the transverse k is part of the result — put TBT.kgrid_Monkhorst_Pack at the very end of input.fdf, and a k convergence test is mandatory
08 Landauer I–VIdenticalApply the same integration formula to the k-averaged T(E)T(E)
09 Finite biasIdenticalTSDE inheritance and the bias-chain file rules carry over unchanged
10 EigenchannelCautionEigenchannels are defined per k\mathbf{k}_\parallel. Analyze at a specific k (usually Γ\Gamma), and when comparing directly with the k-averaged T(E)T(E), state which k the channels belong to
11 Level alignmentCautionA correction premised on a molecular subspace — there is nothing to apply it to in molecule-free bulk and interface systems
12 HPC parallelizationIdenticalThe benefit of k-point and energy-point parallelization becomes even larger

The point of the table is this — reuse the file structure and execution procedure as they are, and re-verify only the k-sampling and its placement rule. Chapters marked "Identical" allow you to copy the input and change only the paths, but for "Caution" chapters, read the corresponding remark first and then proceed.

The recommended learning path also follows from this. Verify the setup with type 1 (a pristine perfect crystal), then move to interfaces and defects (type 2); if a molecular junction is the goal, move on to the metal-electrode molecular junction chapter, where anchoring and scattering-region design are layered on top. Whichever path you take, the principle that results skipping pristine verification cannot be trusted is the same.

Exercises

  1. Au monatomic chain — repeat Chapters 05–07 with the element changed to Au. Confirm the T(EF)=1T(E_F) = 1 plateau in the pristine case, and going below EFE_F, find the energy at which channels of 5d origin open from both the T(E)T(E) staircase and the band structure (Chapter 02 approach), then put the two into correspondence.
  2. k-resolved verification — pick one transverse-periodic system (or merely design it) and plan a plot that overlays the T_k of the snippet above on the energy axis. State as judgment criteria whether each k curve is an integer staircase, whether the weighted average reproduces T_avg, and how T_avg becomes smoother as the number of k increases.
  3. Interface design exercise — write a plan for assembling a metal A/metal B junction device. It must include the transverse cell matching (which side to strain), the initial value of the number of buffer layers per interface and the deliverables for judging convergence, and the role of the two pristine control calculations (A only, B only).

References

  • M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, "Density-functional method for nonequilibrium electron transport", Phys. Rev. B 65, 165401 (2002). DOI 10.1103/PhysRevB.65.165401
  • N. Papior, N. Lorente, T. Frederiksen, A. García, and M. Brandbyge, "Improvements on non-equilibrium and transport Green function techniques: The next-generation TRANSIESTA", Comput. Phys. Commun. 212, 8 (2017). DOI 10.1016/j.cpc.2016.09.022
  • M. P. López Sancho, J. M. López Sancho, and J. Rubio, "Quick iterative scheme for the calculation of transfer matrices", J. Phys. F: Met. Phys. 14, 1205 (1984). DOI 10.1088/0305-4608/14/5/016