SlabTerminator

Enumerate the symmetrically unique slab terminations of a bulk crystal for a given Miller index from crystal symmetry.

Given a bulk pymatgen Structure and a Miller index, SlabTerminator finds every distinct way the crystal can be cleaved along that plane, tells you which terminations are polar vs. nonpolar, and (optionally) builds the ready-to-use slab structures with vacuum.

Why not just use pymatgen?

pymatgen's SlabGenerator.get_slabs() enumerates candidate cleaves by clustering atoms along the normal within a tolerance, builds a full slab for each, and then deduplicates those slabs by comparing them pairwise with StructureMatcher (a tolerance-based lattice-reduction + site-matching comparison). SlabTerminator instead works purely from group theory: it projects the oriented cell's space-group operations onto the surface normal and groups candidate cleaves into symmetry orbits before any vacuum is added (see How it works), so it never builds a slab just to decide uniqueness and never runs a pairwise structure comparison. Two practical consequences:

• It is much faster. On our benchmark of 1352 (structure, Miller) cases, SlabTerminator is roughly 40 to 131× faster than SlabGenerator when building slabs (111× on average) and ~320× faster on average when only counting terminations, because it skips both the per-cleave slab construction and the pairwise StructureMatcher deduplication. In batch mode there is an additional ~26% on top as the bulk symmetry is computed once per material and reused across all its Miller indices.

  

  Per-case speedup (pymatgen SlabGenerator ÷ SlabTerminator) vs. problem size and Miller index, over 1352 cases up to max index 3. Both run at a single fixed layer tolerance (tol = pymatgen ftol = 0.1, no auto scan); matched geometry (≥ 8 Å slab, 10 Å vacuum, max_normal_search=1, centered), best of 3 runs, on an Apple M5 (heat can throttle absolute timings, but methods are timed interleaved so the relative speedup is barely affected) / pymatgen 2026.5.4. "Cleavage planes" = candidate interlayer cleaves before symmetry reduction. Per-case data in scripts/benchmark.csv.

• That speed makes it more robust, by making tolerance sweeps cheap. Both methods share a layer-grouping tolerance (SlabTerminator's tol, pymatgen's ftol) that changes the count: too tight over-splits near-coplanar atoms, too loose merges distinct terminations. Any single tolerance is a guess. Because SlabTerminator's analysis reuses one cached oriented cell and its symmetry operations — no slab rebuilds, no extra spglib calls — sweeping a whole grid of tolerances is nearly free, so it can report the count as a function of tolerance and pick a stable plateau automatically (scan_termination_stability() / tol="auto"; see Choosing the tolerance automatically). Running the same sweep through get_slabs means rebuilding every slab and re-running the pairwise comparison at each tolerance, expensive enough that in practice one picks a single ftol and trusts it.

How it works

The oriented unit cell is periodic along the surface normal. SlabTerminator runs two complementary symmetry analyses, one on each side of adding vacuum:

• Without vacuum → which cleaves are the same slab. The cell's space-group operations are projected onto the 1D coordinate along the surface normal as g → ±g + τ. Candidate interlayer gaps are grouped into orbits under these maps; each orbit is one unique termination. Screw axes, glide planes, and pure c-translations (which only exist while the cell is periodic along the normal) are what relate cleaves recurring at different heights, so this must be done before vacuum is added.
• With vacuum → is a slab polar. Each built slab-with-vacuum is checked for a surviving operation that maps the normal to its negative. If one exists the two faces are equivalent (nonpolar); otherwise the slab is polar. This must be done with vacuum, since the periodic cell can otherwise report a false symmetry through a glide/screw whose translation the vacuum breaks.

Benchmarking

This new method is both faster and more robust than fingerprint-based enumeration (an older version of this code from 2022; unpublished). See scripts/benchmark.py (which stores every per-case count and timing in scripts/benchmark.csv) for a correctness + speed comparison against the old UniqueSlabsGenerator and pymatgen'sSlabGenerator (based on StructureMatcher).

By default the benchmark script recomputes all timings; the --reuse flag pulls whole method groups [new (SlabTerminator), old (old fingerprint), pmg (pymatgen)] from the existing CSV instead of re-measuring them (e.g. --reuse pmg re-times the fast ST methods fresh while keeping the slow, cached pymatgen numbers). Each group keeps its own "measured at" timestamp in the meta sidecar. Using the flag --reuse all produces only the report. scripts/benchmark_analysis.py reads that CSV to plot the speedup-vs-size analysis (scripts/benchmark_speedup_vs_size.png) displayed above and writes scripts/benchmark_summary.md.

A note on AI-assisted development

The original (unpublished) version of this software was written in 2021/22 and used a fingerprint approach based on nearest-neighbor analysis to distinguish terminations. With the help of Claude Code, an entirely new approach was developed from group theory and symmetry operations, and I worked extensively back and forth with AI to ensure its fidelity against both pymatgen and the older fingerprint method, as well as to optimize and analyze the resulting speedup.

Installation

Requires Python ≥ 3.12. The package is on PyPI:

# with uv
uv add slabterminator

# with pip
pip install slabterminator

The only runtime dependency is pymatgen.

To work on the source instead, clone the repo and use uv to install it with its dev dependencies:

git clone https://github.com/d2r2group/slabterminator
cd slabterminator
uv sync

Quick start

from pymatgen.core import Structure
from slabterminator.core import SlabTerminator

structure = Structure.from_file("tests/test-cifs/Fe3C_mp-13154_conventional_standard.cif")

# Analyze the (1, 0, 1) surface.
gen = SlabTerminator(structure, (1, 0, 1))

# Cheap: just enumerate the unique terminations (no slabs built).
for term in gen.get_unique_terminations():
    print(term)
# Termination(gap_index=..., gap_position=0.125, multiplicity=..., symmetric_by_bulk=False)
# ... 4 terminations for Fe3C(101)

# Full: build one slab per unique termination, with vacuum.
result = gen.get_unique_slabs(
    vacuum_size=15.0,          # Angstrom of vacuum along c
    min_slab_thickness=10.0,   # grow the slab until it exceeds this thickness (Angstrom)
    max_normal_search=1,       # search for a more orthogonal output cell
)

print(result.properties.n_unique_terminations)  # 4
print(round(result.properties.surface_area, 2))  # 27.56

for i, entry in enumerate(result.slabs):
    print(i, round(entry.shift, 4),
          "polar" if not entry.is_symmetric_with_vacuum else "nonpolar",
          entry.top_layer_composition, "/", entry.bottom_layer_composition)
    entry.slab.to(filename=f"Fe3C_101_{i}.cif")   # entry.slab is a pymatgen Structure

Output:

0 0.125  polar Fe / Fe
1 0.1844 polar Fe / C
2 0.2292 polar C  / Fe
3 0.4553 polar Fe / Fe

API

SlabTerminator(structure, miller_index, tol=0.1, symprec=0.1, sym_tol=1e-3, slab_symprec=None, tol_scan=None, bulk_symmetry_ops=None)

Constructs the analyzer for one (structure, miller_index) pair. The symmetry analysis runs here, on the cheapest oriented cell, and is independent of how output slabs are later built. Raises ValueError for the (0, 0, 0) index.

tol is the layer c-tolerance (Angstrom) used to group atoms into atomic layers. Pass tol="auto" to have it chosen automatically from a tolerance-stability scan (see Choosing the tolerance automatically below); tol_scan overrides the tolerances swept in that case.

bulk_symmetry_ops is an optional per-material speedup for batch use: pass the bulk's Cartesian space-group operations (the second value from get_sym_distinct_miller_indices_and_symops) and the oriented cell's projected symmetry operations are reconstructed from them in closed form instead of re-running spglib per Miller index — one spglib call per material rather than one per index, for the same result. slabterminator.pipeline wires this through automatically; direct callers can leave it None (the default), which runs spglib on the oriented cell as before.

• get_unique_terminations() → list[Termination], one per unique termination, each with gap_index, gap_position (fractional c of the cleave), multiplicity (number of candidate cleaves that collapsed into it), and symmetric_by_bulk (cheap pre-vacuum face-symmetry estimate). No slab structures are built; this is the fast path.

• get_unique_slabs(...) → UniqueSlabsResult, building one slab per unique termination. Key options:

  • vacuum_size (default 10.0): vacuum thickness in Angstrom.
  • slab_thickness_cells / min_slab_thickness: stack the oriented cell to a fixed number of repeats, or grow it until it exceeds a target thickness in Angstrom.
  • all_unique_terminations_to_top (default False): also emit the flipped counterpart of each polar slab (the other face brought to the top by a true 180° rotation, not a mirror).
  • center_slab (default True): center the slab along c, else leave vacuum on top.
  • max_normal_search (default None): search for a more orthogonal (but thicker) output cell. Affects only slab geometry, not which terminations are found.
  • force_orthogonal_cell (default False): force c orthogonal to the surface plane as a final step (see docstring for caveats).

  The returned UniqueSlabsResult is a NamedTuple:

  • properties: n_unique_terminations, surface_area, and is_surface_symmetric_without_vacuum.
  • settings: the resolved settings actually used (handy for reproducibility), including tol, the layer c-tolerance actually applied (the plateau value when tol="auto"), so a run can be reproduced exactly by passing that float back.
  • oriented_unit_cell: the cell the slabs were built in.
  • slabs: a list of SlabEntry, each with shift, is_symmetric_with_vacuum, top_layer_composition, bottom_layer_composition, slab (a pymatgen Structure), and face ('as_cut' or 'flipped').

regenerate_slabs(shifts, vacuum_size, oriented_unit_cell, ...)

Rebuilds the final slab Structures directly from stored UniqueSlabsResult fields, skipping the whole symmetry analysis; useful for persisting a compact result (shifts + oriented cell + settings) and reconstructing the slabs later.

Choosing the tolerance automatically

The number of unique terminations can depend on the layer c-tolerance tol: too tight over-splits near-coplanar atoms into separate terminations, too loose merges genuinely distinct ones. This is the same shift-enumeration tolerance pymatgen's SlabGenerator exposes as ftol.

scan_termination_stability() sweeps a grid of tolerances and reports how the count varies, reusing the cached oriented cell and symmetry operations (no rebuild, no extra spglib calls, so the whole scan is nearly free):

scan = SlabTerminator(structure, (3, 2, 3)).scan_termination_stability()
print(scan.curve)          # [(0.01, 8), (0.015, 8), ..., (0.1, 3), ...]
print(scan.chosen_tol)     # 0.0612  -- the selected plateau tolerance
print(scan.chosen_count)   # 5
print(scan.is_ambiguous)   # False

Passing tol="auto" to the constructor runs this scan and adopts the selected tolerance, all on the single oriented cell (no second build):

gen = SlabTerminator(structure, (3, 2, 3), tol="auto")
print(gen.tol)                 # 0.0612  (also gen.tol_scan_result -> TolScanResult)
print(gen.get_unique_slabs().settings.tol)   # 0.0612  (recorded for reproducibility)

Selection anchors on the conventional default 0.1 and only overrides it with cause. The count-vs-tol curve is usually a monotone step-down whose two ends are traps (the fine end over-splits, the loose end collapses toward a single termination), so the rule is:

• if 0.1's count is stable (shared with a neighboring tol), keep 0.1 (auto is a no-op for the common case);
• if 0.1 sits on a lone one-tol ledge, pick the widest interior plateau (a run touching neither scan end, excluding both saturations); its tolerance is the geometric mean of the plateau's endpoints;
• if 0.1 is a ledge with no interior plateau, the count is genuinely tolerance-ambiguous: keep 0.1 and set is_ambiguous=True (with a warning).

So tol="auto" never returns a degenerate over-merged or over-split count; where the answer is truly resolution-dependent it says so rather than guessing.

slabterminator.utils

Helpers used by the core, including get_sym_distinct_miller_indices_and_symops(structure, max_index) to enumerate the Miller indices worth analyzing. It returns a (miller_indices, bulk_symmetry_ops) tuple from a single spglib call — the ops let SlabTerminator skip its own per-index spglib call (see bulk_symmetry_ops above); pass them through when looping, or ignore them if you don't need the speedup:

from slabterminator.utils import get_sym_distinct_miller_indices_and_symops

millers, bulk_ops = get_sym_distinct_miller_indices_and_symops(structure, max_index=1)
for miller in millers:
    result = SlabTerminator(structure, miller, bulk_symmetry_ops=bulk_ops).get_unique_slabs()
    print(miller, result.properties.n_unique_terminations)

Batch generation over many materials

SlabTerminator handles one (structure, Miller index) pair. Two higher-level modules build on it for high-throughput datasets: one bulk material in, all its slabs out — and many materials in parallel.

slabterminator.pipeline — one material, all Miller indices

build_slabs_for_material(structure, config) enumerates the symmetrically distinct Miller indices (up to config.max_miller_index), runs SlabTerminator on each, and returns a flat list of records — one per built slab — instead of raising on a bad material (it returns MaterialResult(ok=False, error=...) so a batch can keep going). Parameters are grouped into a single SlabGenConfig rather than a long argument list:

from slabterminator.pipeline import build_slabs_for_material, SlabGenConfig

config = SlabGenConfig(
    max_miller_index=3,
    tol="auto",              # pick each surface's layer tolerance from its plateau
    min_slab_thickness=15.0,
    vacuum_size=15.0,
    center_slab=False,
)
result = build_slabs_for_material(structure, config, material_id="mp-13154")

print(result.ok, result.n_slabs)                 # True 90
for rec in result.records:
    print(rec.miller, rec.term_index, rec.face,
          "polar" if not rec.is_symmetric_with_vacuum else "nonpolar",
          rec.top_layer_composition, "/", rec.bottom_layer_composition)
    # rec.slab is the with-vacuum pymatgen Structure; rec.oriented_unit_cell + rec.shift
    # reconstruct the no-vacuum slab via regenerate_slabs(vacuum_size=0.0).

Each SlabRecord carries the built slab, its shift, oriented_unit_cell, and the per-slab and aggregate properties (is_symmetric_with_vacuum, is_surface_symmetric_without_vacuum, surface_area, n_unique_terminations, the resolved tol and max_normal_search, …). The no-vacuum slab is not stored: the oriented cell plus the shift reconstruct it exactly via regenerate_slabs.

slabterminator.batch — many materials in parallel

run_batch(materials, config, *, on_result, ...) fans build_slabs_for_material out across worker processes and streams each finished MaterialResult to a sink callback. It is scheduler- and output-agnostic: you provide the (id, structure) stream and an on_result writer (CSV, database, in-memory list, …). A material whose worker overruns max_material_seconds (or crashes) is killed and recorded as a failure rather than stalling the run — this is why it uses raw processes rather than a pool, whose futures cannot interrupt a running task.

from slabterminator.batch import run_batch

records = []
def sink(result):
    if result.ok:
        records.extend(result.records)

materials = [("mp-13154", struct_a), ("mp-2657", struct_b)]  # structures or as_dict() forms
n = run_batch(materials, config, n_workers=4,
              max_material_seconds=3 * 3600, on_result=sink)

Passing result_transform= maps each MaterialResult to whatever on_result should receive. For successful materials it runs in the worker process, so heavy per-record serialization (e.g. Structure.as_dict → JSON) is parallelized across workers and only the lightweight payload crosses the results queue — keeping the single parent sink from becoming the bottleneck at high throughput (see scripts/build_slab_dataset_slurm.py, which uses it to stream JSON Lines). Without it, on_result receives the raw MaterialResult.

Passing config_path= writes a self-describing JSON manifest once, before any worker launches, so a dataset records how it was made. It has a versions block (slabterminator, pymatgen, python), a UTC generated_at, a batch block (n_workers, max_material_seconds — the latter determines which materials survive), and config, the resolved SlabGenConfig with library-default None fields filled in.

Testing

uv run pytest

The suite currently runs over 7,700 tests, most of them parametrized across the fixtures below. Fixtures live in tests/test-cifs/: one representative structure for each of the 32 crystallographic point groups (spanning all 7 crystal systems), so the symmetry handling is exercised across the full range of surface symmetries. These are Materials Project conventional standard cells, plus one synthetic polar structure (synthetic_polar_Pna21.cif) for the point group mm2.