Calculating the Schottky barrier

The Schottky barrier is the energy barrier that forms at a metal–semiconductor junction. The barrier height is an important characteristic for charge transport across the junction and a critical parameter in the design of semiconductor contacts, transistors, and rectifying diodes. The p-type Schottky barrier height [math]\displaystyle{ \varphi_{\mathrm{p}} }[/math] describes the barrier seen by holes in the semiconductor valence band, and the n-type Schottky barrier height [math]\displaystyle{ \varphi_{\mathrm{n}} }[/math] describes the barrier seen by electrons in the conduction band. In VASP, the potential alignment methodcitation needed can be used to estimate the Schottky barrier.

Three quantities from bulk calculations, the Fermi energy [math]\displaystyle{ E_{\mathrm{F}} }[/math] of the metal, and the valence-band maximum [math]\displaystyle{ E_{\mathrm{VBM}} }[/math] and band gap [math]\displaystyle{ E_{\mathrm{g}} }[/math] of the semiconductor, are needed for the calculation. As well as the macroscopic electrostatic potential difference [math]\displaystyle{ \Delta\bar{V} }[/math] across the interface, which is extracted from a separate interface slab calculation.

Step-by-step instructions

Important: For the calculation of the p- and n-type Schottky barriers, energies, and potentials from different VASP calculations are compared. It is crucial to set a single value for ENCUT for all calculations.

Step 1: Compute the Fermi energy of the metal.

After ensuring that the bulk geometry is optimized, a separate static calculation should be performed with a dense, Γ-centered k-point mesh and ISMEAR = -5 (thedrahedron method with Blöchl corrections), to get an accurate Fermi energy. It can be taken from the OUTCAR and is also easily accessible via py4vasp:

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

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

Band gaps are notoriously underestimated by standard DFT, so for accurate values, it is usually necessary to use DFT+U, hybrid functionals, or GW calculations. The level of theory to choose depends on the semiconductor of interest. Generally, it is advisable to select BANDGAP = KPOINT for more verbose output, and ideally to make a full bandstructure calculation, since the extrema of the bands will be located at high symmetry lines. The fundamental bandgap [math]\displaystyle{ E_{\mathrm{g}} }[/math] and the valence band maximum [math]\displaystyle{ E_{\mathrm{VBM}} }[/math], can be read from the OUTCAR file, or with py4vasp:

import py4vasp as pv

calc  = pv.Calculation.from_path("./Step2/")

E_g   = calc.bandgap.fundamental()
E_VBM = calc.bandgap.valence_band_maximum()

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

Step 3: Build and relax an interface structure.

Building a lattice-matched metal–semiconductor interface slab is outside the scope of this page. Tools such as the CoherentInterfaceBuilder in pymatgen and the interface construction utilities in ASE can be used to construct the initial geometry. An interface structure for a Schottky barrier calculation should satisfy the following requirements:

• Low lattice mismatch. The supercell should be chosen such that the in-plane periodicities of the two surfaces match to within ~2%, minimising artificial biaxial strain.

• Sufficient slab thickness. Each material should contain enough layers that a bulk-like region exists in the centre.

• No vacuum region. A cell with a vacuum will result in lateral strain, compressing the cell to reduce the surface energy. It is better to have two equivalent interfaces in the cell.

For the relaxation, it is usually prudent to fix the layers in the bulk like regions using selective dynamics in the POSCAR file. Only a couple of layers at the interfaces should be relaxed. Depending on the material pairing and the strain in the interface, it might be beneficial to allow the cell to relax only in the direction normal to the interface plane using the LATTICE_CONSTRAINTS INCAR tag.

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. 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.

The planar average [math]\displaystyle{ \bar{V}(z) }[/math] is obtained by averaging the 3D potential over the in-plane (xy) grid at each [math]\displaystyle{ z }[/math] point. A macroscopic average is then computed by applying a uniform box-car convolution of window length [math]\displaystyle{ L }[/math] equal to the bulk layer period of the materials involved; this removes short-range oscillations within each atomic layer while preserving the long-range potential step across the interface. For Al(111) and Si(111), [math]\displaystyle{ L = 9.35 }[/math] Å covers approximately three Al monolayers or three Si bilayers. Read the potential from vaspout.h5, compute both averages, and plot:

import h5py
import numpy as np
from scipy.ndimage import uniform_filter1d
import matplotlib.pyplot as plt

VASPOUT       = "Step5/vaspout.h5"
WINDOW_LENGTH = 9.35              # macroscopic averaging window (Å)

with h5py.File(VASPOUT, "r") as f:
    hartree = f["/results/potential/hartree"][0]   # shape (NGZ, NGY, NGX)
    ionic   = f["/results/potential/ionic"][0]     # shape (NGZ, NGY, NGX)
    lat     = f["/results/positions/lattice_vectors"][:]

c   = np.linalg.norm(lat[2])         # cell length along interface normal (Å)
NGZ = hartree.shape[0]
dz  = c / NGZ

# planar average: mean over the xy-plane at each z
V_planar = (hartree + ionic).mean(axis=(1, 2))    # shape (NGZ,)  [eV]

# macroscopic average: uniform_filter1d is a box-car convolution of `size` points;
# mode='wrap' enforces periodic boundary conditions along z.
n_window = int(round(WINDOW_LENGTH / dz))
V_macro  = uniform_filter1d(V_planar, size=n_window, mode="wrap")

z = np.arange(NGZ) * dz             # z-coordinates in Å

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(z, V_planar, lw=0.7, color="steelblue", alpha=0.5, label="Planar average")
ax.plot(z, V_macro,  lw=2.0, color="tomato",
        label=f"Macroscopic average (L\u202f=\u202f{WINDOW_LENGTH}\u202fÅ)")
ax.set_xlabel("z (Å)")
ax.set_ylabel("Electrostatic potential (eV)")
ax.set_xlim(0, c)
ax.legend()
plt.tight_layout()
plt.savefig("potential_alignment.png", dpi=150)
plt.show()

Inspect the plot to identify the flat bulk-like plateaus on each side of the interface. The plateau values give [math]\displaystyle{ \bar{V}_{\mathrm{m}} }[/math] and [math]\displaystyle{ \bar{V}_{\mathrm{sc}} }[/math]; adjust the z-range masks below to lie within the flat regions, away from the interface transition zones:

z_frac  = z / c
Al_mask = (z_frac > 0.05) & (z_frac < 0.20)   # bulk-like Al region  (adjust)
Si_mask = (z_frac > 0.42) & (z_frac < 0.58)   # bulk-like Si region  (adjust)

V_Al    = V_macro[Al_mask].mean()
V_Si    = V_macro[Si_mask].mean()
delta_V = V_Al - V_Si

print(f"V̄_Al  = {V_Al:+.4f} eV")
print(f"V̄_Si  = {V_Si:+.4f} eV")
print(f"\u0394V̄     = V̄_Al \u2212 V̄_Si = {delta_V:+.4f} eV")

For the Al(111)/Si(111) interface studied here this gives [math]\displaystyle{ \bar{V}_{\mathrm{m}} = -0.69 }[/math] eV, [math]\displaystyle{ \bar{V}_{\mathrm{sc}} = +1.20 }[/math] eV, and [math]\displaystyle{ \Delta\bar{V} = -1.89 }[/math] eV.

Step 7: Compute the Schottky barrier heights.

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]

# Values from Steps 1, 2, and 6
E_F     =  8.085   # eV  (bulk Al, Step 1b)
E_VBM   =  5.449   # eV  (bulk Si HSE06, Step 2b)
E_g     =  1.160   # eV  (bulk Si HSE06, Step 2b)
delta_V = -1.891   # eV  (interface potential alignment, Step 6)

phi_p = delta_V + E_F - E_VBM
phi_n = E_g - phi_p

print(f"phi_p = {phi_p:.3f} eV   (p-type Schottky barrier)")
print(f"phi_n = {phi_n:.3f} eV   (n-type Schottky barrier)")

For the Al(111)/Si(111) interface this gives [math]\displaystyle{ \varphi_{\mathrm{p}} = 0.745 }[/math] eV and [math]\displaystyle{ \varphi_{\mathrm{n}} = 0.415 }[/math] eV.

Recommendations and advice

(To be added.)

Related tags and articles

(To be added.)