Modeling Catalysis and Surface Adsorption with DFT
Materials Science

Modeling Catalysis and Surface Adsorption with DFT

A practical DFT guide to heterogeneous catalysis: slab models, adsorption energy, the d-band model, Sabatier volcano plots, reaction barriers with NEB, and coverage effects.

Modeling Catalysis and Surface Adsorption with DFT
Photo by Thomas T on Unsplash · View photo

Heterogeneous catalysis — where reactions happen on a solid surface — underpins a vast share of the chemical industry, from ammonia synthesis to automotive exhaust cleanup to green-hydrogen electrolysis. Density functional theory (DFT) has become indispensable for understanding and designing these catalysts, because it can compute the energetics of molecules binding and reacting on surfaces atom by atom. This post lays out the core surface-science workflows DFT enables and the concepts — adsorption energy, the d-band model, volcano plots, and reaction barriers — that turn those numbers into catalyst design principles.

Why DFT Is the Workhorse of Computational Catalysis

Catalysis is fundamentally about energetics on a surface: how strongly reactants bind, how easily bonds break and form, and how readily products release. These are all total-energy differences, exactly what DFT computes. Where experiments see only aggregate turnover rates, DFT resolves individual adsorption sites, intermediates, and transition states.

The payoff is rational catalyst design. Rather than screening alloys and oxides by trial and error, computational groups compute a single descriptor — often an adsorption energy — and use it to rank hundreds of candidate surfaces. This descriptor-based screening, pioneered for electrocatalysis, is one of DFT’s clearest success stories.

Slab Models and Surface Energy

A real surface is a semi-infinite solid terminated by a plane of atoms. DFT represents it with a slab model: a few atomic layers of the crystal, periodic in the surface plane, with a vacuum gap above to separate the slab from its periodic image.

flowchart TB
  subgraph slab [Slab model]
    Vac[Vacuum 12–15 Å]
    L1[Layer 1 — relax]
    L2[Layer 2 — relax]
    L3[Layer 3 — fixed bulk]
    L4[Layer 4 — fixed bulk]
    Vac --- L1
    L1 --- L2
    L2 --- L3
    L3 --- L4
  end
  Ads[Adsorbate on top site / bridge / hollow] --> L1

Key construction choices:

  • Facet / Miller index. Different crystal faces — fcc(111), fcc(100), stepped (211) — expose different atomic arrangements and reactivities. The close-packed (111) face is the common starting point for fcc metals.
  • Slab thickness. Enough layers that the interior recovers bulk-like behavior; typically 4-6 layers, with the bottom layers fixed at bulk positions and the top layers allowed to relax.
  • Vacuum spacing. Usually 12-15 Å, converged so periodic images do not interact; add a dipole correction when adsorbates sit on only one face.

Surface Energy

The surface energy quantifies the cost of creating a surface and controls equilibrium crystal shape (Wulff construction) and which facets dominate. It is computed by comparing the slab energy to the bulk energy per atom:

$$ \gamma = \frac{E_{\mathrm{slab}} - N\, E_{\mathrm{bulk}}}{2A} $$

where \(E_{\mathrm{slab}}\) is the total energy of the relaxed slab, \(E_{\mathrm{bulk}}\) is the bulk energy per atom, \(N\) is the number of atoms in the slab, \(A\) is the surface area of one face, and the factor of 2 accounts for two surfaces (top and bottom) per slab.

Lower surface energy means a more stable, more prevalent facet — information that tells you which surfaces a real catalyst nanoparticle actually exposes.

Adsorption Energy: The Central Quantity

The most important number in surface catalysis is the adsorption energy (or binding energy), the energy released when a molecule or atom binds to the surface:

$$ E_{\mathrm{ads}} = E_{\mathrm{slab+ads}} - E_{\mathrm{slab}} - E_{\mathrm{adsorbate}} $$

A more negative \(E_{\mathrm{ads}}\) means stronger binding. The workflow is to place the adsorbate at candidate sites — top (atop one atom), bridge (between two), hollow (fcc or hcp three-fold sites) — relax each, and identify the most stable configuration. For example, CO on Pt(111) prefers the atop site at low coverage, while H often favors hollow sites.

AdsorbateSurfacePreferred siteBinding trend
HPt(111)fcc hollowmoderate, near-optimal for HER
COPt(111)atop (low coverage)strong, can poison
OPt(111)fcc hollowstrong
HAu(111)weak (Au binds poorly)
NFe(111)high-coordinationstrong (key for ammonia synthesis)

Why does platinum bind hydrogen near-optimally while gold barely binds it at all? The d-band model, developed by Hammer and Norskov, gives a remarkably simple electronic-structure answer that DFT can compute directly.

The idea: when an adsorbate orbital interacts with a transition metal’s d states, it forms bonding and antibonding combinations. How much of the antibonding states sit above the Fermi level (and stay empty) determines binding strength. The key descriptor is the d-band center — the average energy of the metal’s d states relative to the Fermi level.

  • A higher (closer to Fermi) d-band center pushes antibonding states above the Fermi level, leaving them empty, giving stronger adsorbate binding.
  • A lower d-band center fills more antibonding states, weakening binding.

This single quantity, obtained from a projected density of states, explains trends across the transition-metal series and under strain or alloying. Compressive strain lowers the d-band center and weakens binding; tensile strain does the reverse. The d-band model turned catalyst design into a search for surfaces with the right d-band center.

The Sabatier Principle and Volcano Plots

The Sabatier principle states that the best catalyst binds intermediates neither too strongly nor too weakly. Bind too weakly and reactants will not stick or activate; bind too strongly and products poison the surface and will not desorb. The optimum is intermediate.

Plotting catalytic activity against a binding-energy descriptor produces the iconic volcano plot: activity rises as binding strengthens, peaks at the Sabatier optimum, then falls. DFT builds these plots from first principles:

  1. Choose a descriptor (e.g. \(E_{\mathrm{ads}}\) of H* for HER, or O*/OH* for ORR).
  2. Compute that adsorption energy on many candidate surfaces.
  3. Relate activity to the descriptor via a microkinetic or scaling model.
  4. Plot activity vs. descriptor → volcano.
  5. The peak identifies the ideal binding strength and best candidates.

For the hydrogen evolution reaction (HER), the descriptor is the hydrogen adsorption free energy, and the volcano peaks near zero — explaining why platinum, which binds H near-thermoneutrally, is the benchmark catalyst, and motivating cheaper alternatives (MoS2 edges, phosphides) that approach the same binding.

Scaling Relations

A subtlety DFT revealed is that adsorption energies of related intermediates (for example O*, OH*, OOH*) are often linearly correlated — scaling relations. These correlations constrain how good any single catalyst can be and define the theoretical activity limits that catalyst designers try to break by disrupting the scaling.

Reaction Barriers with the Nudged Elastic Band Method

Thermodynamics (adsorption energies) tells you what is stable; kinetics tells you how fast. Reaction rates depend on activation barriers at transition states, and the standard DFT tool to find them is the nudged elastic band (NEB) method, especially its climbing-image variant.

NEB connects a known initial state (reactants adsorbed) and final state (products adsorbed) with a chain of intermediate images, then relaxes them to trace the minimum-energy path over the potential energy surface. The highest point is the transition state; its energy relative to the initial state is the barrier.

$$ E_a = E(\mathrm{TS}) - E(\mathrm{initial}) $$

A classic application is CO oxidation (CO* + O* → CO2) on platinum-group metals, where NEB barriers explain the light-off temperature of automotive catalytic converters. Barriers feed directly into rate constants via transition-state theory, closing the loop from electronic structure to measurable kinetics. Quantum ESPRESSO ships an NEB engine for exactly these calculations; see the Quantum ESPRESSO documentation.

Coverage Effects and Realistic Conditions

Textbook adsorption energies assume a single isolated molecule on a large surface — the low-coverage limit. Real catalysts operate at high coverage where adsorbates crowd and interact.

  • Lateral interactions. Neighboring adsorbates repel (or occasionally attract), so binding weakens as coverage increases. Computing \(E_{\mathrm{ads}}\) at 1/9, 1/4, 1/2, and full monolayer coverage maps this out.
  • Adsorbate ordering. Dense layers form ordered superstructures (for example the (2x2) and c(4x2) CO overlayers on Pt(111)) that DFT reproduces in appropriately sized supercells.
  • Operating conditions. Coupling DFT energies to a thermodynamic reservoir via ab initio thermodynamics predicts the equilibrium coverage as a function of temperature and partial pressure, bridging the gap between 0 K calculations and reactor conditions.

Capturing coverage is where supercell size matters most: resolving dilute coverages or complex overlayers on stepped surfaces can require hundreds to thousands of atoms.

From Descriptors to Design: The Full Workflow

A complete computational catalysis study weaves these pieces together:

  1. Build the surface — choose facet, converge slab thickness and vacuum, compute surface energy.
  2. Adsorption — find stable sites and adsorption energies for reactants, intermediates, products.
  3. Electronic analysis — extract the d-band center and rationalize binding trends.
  4. Screening — build volcano plots and apply scaling relations across candidates.
  5. Kinetics — run NEB for the rate-limiting steps.
  6. Realism — add coverage effects and ab initio thermodynamics.

For more on the underlying methods and engines, see our computational engines page and related posts on the blog.

Run it on Simatra

Surface science is computationally intensive: slab convergence, site searches, NEB paths with many images, and high-coverage supercells all multiply the number of expensive relaxations a project needs. Simatra runs these DFT workflows on GPU-accelerated clusters built around the GPU-Opt-V2 instance, delivering up to 5x faster convergence and supporting supercells up to roughly 2,000 atoms — enough for stepped surfaces, alloy facets, and dense adsorbate overlayers. Run the standard Quantum ESPRESSO stack or our native KRONOS engine, and scale from a single adsorption energy to a full volcano-plot screen. Start with a free trial and $100 in credits at app.simatra.io.