Calculating the Schottky barrier

The Schottky barrier is the energy barrier that forms at a metal–semiconductor junction when the two materials are brought into contact and their Fermi levels align at equilibrium. Its height governs charge transport across the junction and is a critical parameter in the design of semiconductor contacts, transistors, and rectifying diodes. Two complementary quantities characterise the barrier: the p-type Schottky barrier height [math]\displaystyle{ \varphi_{\mathrm{p}} }[/math], which describes the barrier seen by holes in the semiconductor valence band, and the n-type Schottky barrier height [math]\displaystyle{ \varphi_{\mathrm{n}} }[/math], which describes the barrier seen by electrons in the conduction band. This page describes how to compute both quantities for the Al(111)/Si(111) interface using the potential alignment methodcitation needed with VASP and py4vasp.

The potential alignment method obtains the Schottky barrier heights by combining three quantities from bulk calculations — the metal Fermi energy [math]\displaystyle{ E_{\mathrm{F}} }[/math], the semiconductor valence-band maximum [math]\displaystyle{ E_{\mathrm{VBM}} }[/math], and the semiconductor band gap [math]\displaystyle{ E_{\mathrm{g}} }[/math] — with the macroscopic electrostatic potential difference [math]\displaystyle{ \Delta\bar{V} }[/math] across the interface, which is extracted from a separate interface slab calculation. In VASP, the [math]\displaystyle{ \mathbf{G}=0 }[/math] Fourier component of the Hartree potential is set to zero by convention, so all single-particle eigenvalues are already referenced to the average electrostatic potential of their respective unit cell. As a consequence, [math]\displaystyle{ E_{\mathrm{F}} }[/math] and [math]\displaystyle{ E_{\mathrm{VBM}} }[/math] can be read directly from the bulk OUTCAR files without any further potential referencing, and only the interface calculation requires the electrostatic potential written to LOCPOT.

Step-by-step instructions

Step 1: Compute the Fermi energy of the metal using bulk Al. Steps 1 and 2 are fully independent and can be submitted simultaneously.

Step 1a: Geometry optimization of bulk Al

Perform a full geometry optimization of the bulk metal unit cell:

IBRION = 1 ; ISIF = 7 ; NSW = 50
ISMEAR = 1 ; SIGMA = 0.1
ENCUT = 300 ; EDIFF = 1E-6 ; EDIFFG = -0.01

ISIF = 7 allows the cell shape, volume, and ionic positions to relax simultaneously. ISMEAR = 1 (Methfessel-Paxton) is the correct smearing scheme for metals.

Mind: Use the same ENCUT in all calculations of this workflow to ensure a consistent potential alignment.

Step 1b: Static calculation of bulk Al

Copy CONTCAR to POSCAR and run a single-point calculation:

ISMEAR = -5
ENCUT = 300 ; PREC = Accurate ; EDIFF = 1E-8

ISMEAR = -5 (tetrahedron method with Blöchl corrections) gives an accurate Fermi energy from a static calculation on a regular k-mesh. Extract [math]\displaystyle{ E_{\mathrm{F}} }[/math] with py4vasp:

import py4vasp as pv
calc = pv.Calculation.from_path("Al_bulk_static")
E_F = calc.dos.read()["fermi_energy"]
print(f"E_F = {E_F:.4f} eV")

For bulk Al this gives [math]\displaystyle{ E_{\mathrm{F}} = 8.085 }[/math] eV.

Step 2: Compute the band gap and valence-band maximum of the semiconductor using bulk Si.

Step 2a: Geometry optimization of bulk Si

Use the PBEsol functional for an accurate equilibrium lattice constant:

GGA = PS
IBRION = 1 ; ISIF = 7 ; NSW = 50
ISMEAR = -5 ; SIGMA = 0.05
ENCUT = 300 ; EDIFF = 1E-6 ; EDIFFG = -0.01

Step 2b: Band-structure calculation of bulk Si with HSE06

The band gap and valence-band maximum are computed with the HSE06 hybrid functional, which gives accurate band gaps for semiconductors. This requires a hybrid band-structure calculation; follow Band-structure calculation using hybrid functionals to set up the KPOINTS and KPOINTS_OPT files. The key INCAR settings are:

LHFCALC = .TRUE. ; HFSCREEN = 0.2
GGA = PE ; HFRCUT = -1
ALGO = Normal ; ISMEAR = -5
ENCUT = 300 ; PREC = Accurate ; EDIFF = 1E-6

GGA = PE (PBE) is the correct base functional for HSE06, not PBEsol. HFRCUT = -1 avoids discontinuities in the band structure for gapped systems. Extract [math]\displaystyle{ E_{\mathrm{g}} }[/math] and [math]\displaystyle{ E_{\mathrm{VBM}} }[/math] with py4vasp:

import py4vasp as pv
import numpy as np

calc  = pv.Calculation.from_path("Si_bulk_HSE_bands")
bg    = calc.bandgap.read()
band  = calc.band.read()

E_g   = bg["fundamental"]
E_F   = bg["fermi_energy"]
E_VBM = (band["bands"][0] + E_F)[band["occupations"][0] > 0.5].max()

print(f"E_g   = {E_g:.4f} eV")
print(f"E_VBM = {E_VBM:.4f} eV")

For bulk Si this gives [math]\displaystyle{ E_{\mathrm{g}} = 1.160 }[/math] eV and [math]\displaystyle{ E_{\mathrm{VBM}} = 5.449 }[/math] eV.

Step 3: Build the Al(111)/Si(111) interface structure.

Building a lattice-matched metal–semiconductor interface slab is outside the scope of this page. Tools such as the CoherentInterfaceBuilder in pymatgencitation needed and the interface construction utilities in ASEcitation needed can be used to construct the initial geometry. The structure must satisfy the following requirements before it is passed to Step 4:

• Low lattice mismatch. A coincidence-site lattice supercell should be chosen such that the in-plane periodicities of the two surfaces match to within ~2%, minimising artificial biaxial strain. For Al(111) and Si(111), a 4×4 Al supercell matched to a 3×3 Si supercell achieves a mismatch of 0.9%.

• Sufficient slab thickness. Each material should contain enough layers that a bulk-like region, free of interface perturbations, exists in the centre. A minimum of 8–10 monolayers per material is a practical starting point; convergence with slab thickness should be verified. The slab used here contains 13 Al(111) monolayers and 12 Si(111) monolayers (6 bilayers).

• Selective dynamics. Only the layers immediately adjacent to the interface (2–3 layers on each side) should be free to relax in Step 4; the remaining inner layers are frozen. This substantially reduces the cost of the relaxation.

• Reasonable initial interface geometry. The initial separation between the last metal layer and the first semiconductor layer should be set to a physically plausible value (2–3 Å, close to the equilibrium Al–Si bond length) to ensure smooth convergence.

A pre-optimized POSCAR for this interface will be provided here; readers starting from that file may proceed directly to Step 5.

Step 4: Relax the interface.

The relaxation is performed in two passes. In the first pass, only the ionic positions are updated with the cell held fixed (ISIF = 2). After convergence, copy CONTCAR to POSCAR and run a second pass with full cell relaxation (ISIF = 3), which relieves any residual stress from the lattice mismatch. LREAL = A is recommended throughout for this large supercell.

Step 4a: Ionic relaxation (fixed cell)

ALGO = Fast ; LREAL = A
ISMEAR = 0 ; SIGMA = 0.1
ENCUT = 300 ; PREC = Normal
EDIFF = 1E-6 ; EDIFFG = -0.01
IBRION = 2 ; ISIF = 2 ; NSW = 200

ISMEAR = 0 (Gaussian smearing) is the safe choice for the interface supercell, which contains both metallic and semiconducting regions. LREAL = A uses real-space projection for the non-local operations, which is significantly faster for cells with more than ~100 atoms.

Step 4b: Full cell relaxation

Copy CONTCAR to POSCAR and re-run with cell degrees of freedom enabled:

IBRION = 1 ; ISIF = 3

All other tags are unchanged from Step 4a.

Step 5: Compute the electrostatic potential of the interface.

To extract the potential alignment, a static calculation is performed on the relaxed interface with WRT_POTENTIAL = hartree ionic, which writes the ionic and Hartree contributions to the electrostatic potential into vaspout.h5 for direct access via py4vasp. The electrostatic (ionic + Hartree) part of the potential is the correct choice for potential alignment; the exchange-correlation contribution present in the total potential (LVTOT) must be excluded.

ALGO = Fast ; LREAL = A
ISMEAR = -5
ENCUT = 300 ; PREC = Accurate ; EDIFF = 1E-8
WRT_POTENTIAL = hartree ionic

Step 6: Extract [math]\displaystyle{ \Delta\bar{V} }[/math] from the electrostatic potential.

(py4vasp script using calculation.potential for planar and macroscopic averaging along the interface normal, figure of the macroscopic average potential, and procedure for reading [math]\displaystyle{ \Delta\bar{V} }[/math] from the flat bulk-like regions on each side to be added.)

Step 7: Compute the Schottky barrier height.

With [math]\displaystyle{ E_{\mathrm{F}} }[/math], [math]\displaystyle{ E_{\mathrm{VBM}} }[/math], [math]\displaystyle{ E_{\mathrm{g}} }[/math], and [math]\displaystyle{ \Delta\bar{V} }[/math] in hand, the potential alignment difference [math]\displaystyle{ \Delta\bar{V} = \bar{V}_{\mathrm{m}} - \bar{V}_{\mathrm{sc}} }[/math] (positive if the metal side has the higher average potential) connects the bulk reference frames. The p-type and n-type Schottky barrier heights are then

[math]\displaystyle{ \varphi_{\mathrm{p}} = \Delta\bar{V} + E_{\mathrm{F}} - E_{\mathrm{VBM}}, }[/math]

[math]\displaystyle{ \varphi_{\mathrm{n}} = E_{\mathrm{g}} - \varphi_{\mathrm{p}}. }[/math]

(py4vasp script that loads all three calculations, assembles the four input quantities, and evaluates the two equations above to be added. Worked numerical example with the Al/Si values to be added.)

Recommendations and advice

(To be added.)

Related tags and articles

(To be added.)