diff --git a/docs/source/user_guide/benchmarks/surfaces.rst b/docs/source/user_guide/benchmarks/surfaces.rst
index 780118c7c..3cbf32110 100644
--- a/docs/source/user_guide/benchmarks/surfaces.rst
+++ b/docs/source/user_guide/benchmarks/surfaces.rst
@@ -145,3 +145,56 @@ Reference data:
   * S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. Chevrier, K. A. Persson, G. Ceder, "Python Materials Genomics (pymatgen): A Robust, Open-Source Python Library for Materials Analysis," Comput. Mater. Sci., 2013, 68, 314–319. https://doi.org/10.1016/j.commatsci.2012.10.028
 
 * Tran et al. relaxed the slabs using spin-polarized PBE calculations performed in VASP, with a cutoff energy of 400 eV.
+
+Cleavage Energy
+===============
+
+Summary
+-------
+
+Performance in predicting cleavage energies for 36,718 surface configurations
+across a wide range of materials and Miller indices.
+
+Metrics
+-------
+
+1. Cleavage energy MAE
+
+Accuracy of cleavage energy predictions compared to DFT reference values.
+
+For each surface, the cleavage energy is calculated as
+``(E_slab - thickness_ratio * E_bulk) / (2 * A)``, where ``E_slab`` and
+``E_bulk`` are single-point energies of the slab and the lattice-matched bulk
+unit cell, ``thickness_ratio`` is the number of bulk unit cells in the slab
+thickness, and ``A`` is the surface area. Results are reported in meV/A^2.
+The mean absolute error is computed over all 36,718 surfaces.
+
+2. Cleavage energy RMSE
+
+Root mean squared error of cleavage energy predictions across all surfaces.
+
+Computational cost
+------------------
+
+Medium: benchmark involves only single-point calculations, but for 36,718 slab-bulk pairs. Takes roughly 5-20 minutes on GPU or a few hours on CPUs.
+
+Data availability
+-----------------
+
+Input data:
+
+* Surface configurations were obtained from the Materials Project, covering
+  3,699 unique bulk materials with multiple Miller indices and terminations
+  per material. The original unfiltered data source is available at
+  Zenodo (DOI: 10.5281/zenodo.10381505).
+
+Reference data:
+
+* DFT cleavage energies calculated using PBE functional.
+
+Publication:
+
+* A. Mehdizadeh and P. Schindler, "Surface stability modeling with universal
+  machine learning interatomic potentials: a comprehensive cleavage energy
+  benchmarking study," Mach. Learn.: Sci. Technol., 2025.
+  https://iopscience.iop.org/article/10.1088/3050-287X/ae1408
diff --git a/ml_peg/analysis/surfaces/cleavage_energy/analyse_cleavage_energy.py b/ml_peg/analysis/surfaces/cleavage_energy/analyse_cleavage_energy.py
new file mode 100644
index 000000000..af5d14b8a
--- /dev/null
+++ b/ml_peg/analysis/surfaces/cleavage_energy/analyse_cleavage_energy.py
@@ -0,0 +1,225 @@
+"""Analyse cleavage energy benchmark."""
+
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from ml_peg.analysis.utils.decorators import build_table, plot_density_scatter
+from ml_peg.analysis.utils.utils import (
+    load_metrics_config,
+    mae,
+    write_density_trajectories,
+)
+from ml_peg.app import APP_ROOT
+from ml_peg.calcs import CALCS_ROOT
+from ml_peg.models.get_models import get_model_names
+from ml_peg.models.models import current_models
+
+MODELS = get_model_names(current_models)
+CALC_PATH = CALCS_ROOT / "surfaces" / "cleavage_energy" / "outputs"
+OUT_PATH = APP_ROOT / "data" / "surfaces" / "cleavage_energy"
+
+METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml")
+DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config(
+    METRICS_CONFIG_PATH
+)
+
+EV_TO_MEV = 1000.0
+
+
+def compute_cleavage_energy(
+    slab_energy: float,
+    bulk_energy: float,
+    thickness_ratio: float,
+    area_slab: float,
+) -> float:
+    """
+    Compute cleavage energy from slab and bulk energies.
+
+    Parameters
+    ----------
+    slab_energy
+        Total energy of the slab.
+    bulk_energy
+        Total energy of the lattice-matched bulk unit cell.
+    thickness_ratio
+        Number of bulk unit cells in the slab thickness.
+    area_slab
+        Surface area of the slab in Angstrom^2.
+
+    Returns
+    -------
+    float
+        Cleavage energy in eV/Angstrom^2.
+    """
+    return (slab_energy - thickness_ratio * bulk_energy) / (2 * area_slab)
+
+
+@pytest.fixture
+def cleavage_energies() -> dict[str, dict[str, list]]:
+    """
+    Get cleavage energies for all systems in meV/A^2.
+
+    Also saves per-system errors for the distribution plot.
+
+    Returns
+    -------
+    dict[str, dict[str, list]]
+        Dictionary of model names to ``{"ref": [...], "pred": [...]}`` in meV/A^2.
+    """
+    results = {mlip: {"ref": [], "pred": []} for mlip in MODELS}
+    ref_stored = False
+    stored_ref = []
+
+    for model_name in MODELS:
+        model_dir = CALC_PATH / model_name
+        if not model_dir.exists():
+            continue
+
+        model_pred = []
+        model_ref = []
+
+        for xyz_file in sorted(model_dir.glob("*.xyz"), key=lambda p: int(p.stem)):
+            slab = read(xyz_file)
+
+            pred_ce = (
+                compute_cleavage_energy(
+                    slab.info["slab_energy"],
+                    slab.info["bulk_energy"],
+                    slab.info["thickness_ratio"],
+                    slab.info["area_slab"],
+                )
+                * EV_TO_MEV
+            )
+            model_pred.append(pred_ce)
+
+            if not ref_stored:
+                model_ref.append(slab.info["ref_cleavage_energy"] * EV_TO_MEV)
+
+        if model_pred:
+            results[model_name]["pred"] = model_pred
+            if not ref_stored:
+                stored_ref = model_ref
+            results[model_name]["ref"] = stored_ref
+            ref_stored = True
+
+    return results
+
+@pytest.fixture
+@plot_density_scatter(
+    filename=OUT_PATH / "figure_cleavage_energies.json",
+    title="Cleavage Energies",
+    x_label="Predicted cleavage energy / meV/Ų",
+    y_label="Reference cleavage energy / meV/Ų",
+)
+def cleavage_density(
+    cleavage_energies: dict[str, dict[str, list]],
+) -> dict[str, dict]:
+    """
+    Build density scatter inputs and write density trajectories.
+
+    Parameters
+    ----------
+    cleavage_energies
+        Reference and predicted cleavage energies per model.
+
+    Returns
+    -------
+    dict[str, dict]
+        Mapping of model names to density-plot payloads.
+    """
+    label_list = [
+        f.stem
+        for f in sorted(
+            (CALC_PATH / MODELS[0]).glob("*.xyz"), key=lambda p: int(p.stem)
+        )
+    ]
+
+    density_inputs: dict[str, dict] = {}
+    for model_name in MODELS:
+        preds = cleavage_energies[model_name]["pred"]
+        refs = cleavage_energies[model_name]["ref"]
+        density_inputs[model_name] = {
+            "ref": refs,
+            "pred": preds,
+            "meta": {"system_count": len([v for v in preds if v is not None])},
+        }
+        if preds:
+            write_density_trajectories(
+                labels_list=label_list,
+                ref_vals=refs,
+                pred_vals=preds,
+                struct_dir=CALC_PATH / model_name,
+                traj_dir=OUT_PATH / model_name / "density_traj",
+                struct_filename_builder=lambda label: f"{label}.xyz",
+            )
+    return density_inputs
+@pytest.fixture
+def cleavage_mae(cleavage_energies: dict[str, list]) -> dict[str, float]:
+    """
+    Get mean absolute error for cleavage energies.
+
+    Parameters
+    ----------
+    cleavage_energies
+        Dictionary of reference and predicted cleavage energies.
+
+    Returns
+    -------
+    dict[str, float]
+        MAE for each model.
+    """
+    results = {}
+    for model_name in MODELS:
+        ref_vals = cleavage_energies[model_name]["ref"]
+        pred_vals = cleavage_energies[model_name]["pred"]
+        if pred_vals:
+            results[model_name] = mae(ref_vals, pred_vals)
+        else:
+            results[model_name] = None
+    return results
+
+
+@pytest.fixture
+@build_table(
+    filename=OUT_PATH / "cleavage_energy_metrics_table.json",
+    metric_tooltips=DEFAULT_TOOLTIPS,
+    thresholds=DEFAULT_THRESHOLDS,
+    weights=DEFAULT_WEIGHTS,
+)
+def metrics(cleavage_mae: dict[str, float]) -> dict[str, dict]:
+    """
+    Get all cleavage energy metrics.
+
+    Parameters
+    ----------
+    cleavage_mae
+        Mean absolute errors for all models.
+
+    Returns
+    -------
+    dict[str, dict]
+        Metric names and values for all models.
+    """
+    return {"MAE": cleavage_mae}
+
+
+def test_cleavage_energy(
+    metrics: dict[str, dict],
+    cleavage_density: dict[str, dict],
+) -> None:
+    """
+    Run cleavage energy analysis.
+
+    Parameters
+    ----------
+    metrics
+        All cleavage energy metrics.
+    cleavage_density
+        Density-scatter inputs for all models (drives saved plots).
+    """
+    return
diff --git a/ml_peg/analysis/surfaces/cleavage_energy/metrics.yml b/ml_peg/analysis/surfaces/cleavage_energy/metrics.yml
new file mode 100644
index 000000000..c5839cc62
--- /dev/null
+++ b/ml_peg/analysis/surfaces/cleavage_energy/metrics.yml
@@ -0,0 +1,7 @@
+metrics:
+  MAE:
+    good: 0.0
+    bad: 10.0
+    unit: meV/A^2
+    tooltip: Mean Absolute Error of cleavage energies
+    level_of_theory: PBE
diff --git a/ml_peg/app/surfaces/cleavage_energy/app_cleavage_energy.py b/ml_peg/app/surfaces/cleavage_energy/app_cleavage_energy.py
new file mode 100644
index 000000000..b705c5762
--- /dev/null
+++ b/ml_peg/app/surfaces/cleavage_energy/app_cleavage_energy.py
@@ -0,0 +1,95 @@
+"""Run cleavage energy app."""
+
+from __future__ import annotations
+
+from dash import Dash
+from dash.html import Div
+
+from ml_peg.app import APP_ROOT
+from ml_peg.app.base_app import BaseApp
+from ml_peg.app.utils.build_callbacks import plot_from_table_cell, struct_from_scatter
+from ml_peg.app.utils.load import collect_traj_assets, read_density_plot_for_model
+from ml_peg.models.get_models import get_model_names
+from ml_peg.models.models import current_models
+
+MODELS = get_model_names(current_models)
+BENCHMARK_NAME = "Cleavage Energy"
+DOCS_URL = (
+    "https://ddmms.github.io/ml-peg/user_guide/benchmarks/surfaces.html#cleavage-energy"
+)
+DATA_PATH = APP_ROOT / "data" / "surfaces" / "cleavage_energy"
+
+
+class CleavageEnergyApp(BaseApp):
+    """Cleavage energy benchmark app layout and callbacks."""
+
+    def register_callbacks(self) -> None:
+        """Register callbacks to app."""
+        density_plots: dict[str, dict] = {}
+        for model in MODELS:
+            density_graph = read_density_plot_for_model(
+                filename=DATA_PATH / "figure_cleavage_energies.json",
+                model=model,
+                id=f"{BENCHMARK_NAME}-{model}-density",
+            )
+            if density_graph is not None:
+                density_plots[model] = {"MAE": density_graph}
+
+        plot_from_table_cell(
+            table_id=self.table_id,
+            plot_id=f"{BENCHMARK_NAME}-figure-placeholder",
+            cell_to_plot=density_plots,
+        )
+
+        struct_trajs = collect_traj_assets(
+            data_path=DATA_PATH,
+            assets_prefix="assets/surfaces/cleavage_energy",
+            models=MODELS,
+            traj_dirname="density_traj",
+            suffix=".extxyz",
+        )
+        for model in struct_trajs:
+            struct_from_scatter(
+                scatter_id=f"{BENCHMARK_NAME}-{model}-density",
+                struct_id=f"{BENCHMARK_NAME}-struct-placeholder",
+                structs=struct_trajs[model],
+                mode="traj",
+            )
+
+
+def get_app() -> CleavageEnergyApp:
+    """
+    Get cleavage energy benchmark app layout and callback registration.
+
+    Returns
+    -------
+    CleavageEnergyApp
+        Benchmark layout and callback registration.
+    """
+
+    return CleavageEnergyApp(
+        name=BENCHMARK_NAME,
+        description=(
+            "Performance in predicting cleavage energies for "
+            "36,718 surface configurations."
+        ),
+        docs_url=DOCS_URL,
+        table_path=DATA_PATH / "cleavage_energy_metrics_table.json",
+        extra_components=[
+            Div(id=f"{BENCHMARK_NAME}-figure-placeholder"),
+            Div(id=f"{BENCHMARK_NAME}-struct-placeholder"),
+        ],
+    )
+
+
+if __name__ == "__main__":
+    full_app = Dash(
+        __name__,
+        assets_folder=DATA_PATH.parent.parent,
+    )
+
+    cleavage_energy_app = get_app()
+    full_app.layout = cleavage_energy_app.layout
+    cleavage_energy_app.register_callbacks()
+
+    full_app.run(port=8056, debug=True)
diff --git a/ml_peg/calcs/surfaces/cleavage_energy/calc_cleavage_energy.py b/ml_peg/calcs/surfaces/cleavage_energy/calc_cleavage_energy.py
new file mode 100644
index 000000000..fa8507e1f
--- /dev/null
+++ b/ml_peg/calcs/surfaces/cleavage_energy/calc_cleavage_energy.py
@@ -0,0 +1,77 @@
+"""Run calculations for cleavage energy benchmark."""
+
+from __future__ import annotations
+
+from copy import copy
+import json
+from pathlib import Path
+from typing import Any
+
+from ase.io import read, write
+from tqdm import tqdm
+import pytest
+
+from ml_peg.calcs.utils.utils import download_s3_data
+from ml_peg.models.get_models import load_models
+from ml_peg.models.models import current_models
+
+MODELS = load_models(current_models)
+
+OUT_PATH = Path(__file__).parent / "outputs"
+
+
+@pytest.mark.parametrize("mlip", MODELS.items())
+def test_cleavage_energy(mlip: tuple[str, Any]) -> None:
+    """
+    Run cleavage energy benchmark calculations.
+
+    For each surface configuration, single-point energies are computed for
+    the slab and the lattice-matched bulk. Results are saved as a single
+    JSON file per model containing only energies and system identifiers,
+    avoiding redundant storage of atomic coordinates.
+
+    Parameters
+    ----------
+    mlip
+        Name of model and model to get calculator.
+    """
+    model_name, model = mlip
+    calc = model.get_calculator()
+
+    data_dir = (
+        download_s3_data(
+            key="inputs/surfaces/cleavage_energy/cleavage_energy.zip",
+            filename="cleavage_energy.zip",
+        )
+        / "cleavage_energy"
+    )
+
+    write_dir = OUT_PATH / model_name
+    write_dir.mkdir(parents=True, exist_ok=True)
+
+    idx = 0
+    for mpid_dir in tqdm(sorted(d for d in data_dir.iterdir() if d.is_dir())):
+        for xyz_file in sorted(mpid_dir.glob("*.xyz")):
+            structs = read(xyz_file, index=":")
+            slab, bulk = structs[0], structs[1]
+
+            slab.calc = copy(calc)
+            slab_energy = float(slab.get_potential_energy())
+
+            bulk.calc = copy(calc)
+            bulk_energy = float(bulk.get_potential_energy())
+
+            slab.info.update(
+                {
+                    "slab_energy": slab_energy,
+                    "bulk_energy": bulk_energy,
+                    "area_slab": float(slab.info["area_slab"]),
+                    "thickness_ratio": float(slab.info["thickness_ratio"]),
+                    "ref_cleavage_energy": float(slab.info["ref_cleavage_energy"]),
+                    "mpid": slab.info["mpid"],
+                    "miller": slab.info["miller"],
+                    "term": int(slab.info["term"]),
+                }
+            )
+            write(write_dir / f"{idx}.xyz", slab, format="extxyz")
+            idx += 1