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. For plain slab generation (symmetrize=False; left plot), on our benchmark of 30,302 (material, Miller-index) cases across 1000 materials, SlabTerminator is, per case, roughly 15× (index 1) to 92× (index 3) faster than SlabGenerator when building slabs, a geometric-mean ~62× across the full dataset (median ~75×; the tail-weighted whole-workload aggregate is ~97×). And ~108× faster when only counting terminations, because it skips both the per-cleave slab construction and the pairwise StructureMatcher deduplication. All per-case figures here are geometric means of the per-case ratios. In batch mode there is a further ~1.7× on top (up to 2.83×) as the bulk symmetry is computed once per material and reused across all its Miller indices.

  When the slabs are being symmetrized (symmetrize=True), the per-case generation speedup over pymatgen is even more pronounced: geometric-mean roughly 19× (index 1) to 197× (index 3) per case (~113× across the full dataset; the tail-weighted aggregate is ~274×). This per-case figure is higher than the symmetrize=False speedup because pymatgen's atom-removal symmetrizer has pathological worst cases (its slowest polar level-3 orientations run into seconds apiece), whereas SlabTerminator trims straight to a verified crystal flip center in closed form. pymatgen additionally trims each polar slab atom-by-atom and re-runs its symmetry analysis (detailed speedup breakdown, see right plot).

  Both benchmarks are run with scripts/benchmark.py and the full list of mpids used for benchmarking is stored in benchmark_HPC_<date>.meta.md.

  

  Per-case speedup (pymatgen SlabGenerator ÷ SlabTerminator) vs. problem size and Miller index, up to max index 3. Left: plain generation (symmetrize=False), right: symmetrized generation (symmetrize=True); all 30,302 cases completed for both tools (300 s timeout, none excluded). The top scatter shows every case; the two bar panels show the geometric-mean speedup per bin with multiplicative geometric-standard-deviation whiskers. The two columns share matched y-axes per panel, so the higher symmetrize=True speedups are directly readable against the left. Both compare against SlabTerminator 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 (best of 2 for pymatgen), on an Intel Xeon Gold 6448Y HPC node (all methods timed interleaved per case, so the relative speedup is unaffected by any drift) / pymatgen 2026.5.4. "Cleavage planes" = candidate interlayer cleaves before symmetry reduction. Per-case data in date-stamped scripts/benchmark_<date>.csv and scripts/benchmark_symmetrize_<date>.csv.

  This holds up end-to-end, too. The same 1000-material benchmark set, generated to a Miller index of 3 in full batch mode (10 workers on Northeastern University's Explorer HPC Cluster, build + serialization + JSONL I/O), took 67.5 s with SlabTerminator versus 10121.4 s with pymatgen for symmetrize=False (72k slabs each), a ~150× whole-batch speedup on real production hardware at a matched fixed layer tolerance (tol = ftol = 0.1). With symmetrize=True the same build took 102.8 s versus 11446.1 s (~111×), even though SlabTerminator emits (and serializes to disk) ~1.66× more slabs than pymatgen there (143k vs. 86k), because it keeps the polar parents and unsymmetrizable Tasker-III terminations pymatgen drops (SlabTerminator can reproduce pymatgen's behavior of dropping every polar slab by passing symmetrize="drop_all_polar" instead of symmetrize=True, keeping only the slabs it can make symmetric). The two SLURM entry points used are scripts/build_slab_dataset_slurm.py and scripts/build_slab_dataset_slurm_pymatgen.py, matched on everything but the generator (run-wide details in scripts/benchmark_HPC_<date>.meta.md).

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

• Its polarity verdict is cell-orientation-invariant. To decide whether a slab is nonpolar, SlabTerminator applies the bulk flipping operations directly to the built slab's atoms and tests self-coincidence in the surface {a, b} lattice, so the answer depends only on the surface, not on the choice of max_normal_search. Running the full SpacegroupAnalyzer on the slab (pymatgen's route) instead misses a genuine flip when that choice leaves the cell's c vector oblique, making the verdict max_normal_search- dependent: e.g. 4mm(1,1,0), 32(1,0,0), 32(1,0,1), 3m(1,1,0), -6m2(1,1,0) and Fe3C(1,0,2) all read polar at the cheap cell but nonpolar once it is orthogonalized. Because that verdict is also thickness-invariant, SlabTerminator checks it on the minimal single-cell oriented cell rather than the full slab, so for a 15 Å-thick slab it is ~5× faster per termination than the full-slab SpacegroupAnalyzer call pymatgen must run. (spglib remains available as a cross-check oracle via slab_symmetry_method="spglib".)

• A full surface-polarity classification is built in. Because the orbits already come from the crystal's projected symmetry operations, that same machinery classifies each orientation's polarity at no extra cost: whether a flipping operation exists at all (has_flip_op_by_bulk, equivalently whether any nonpolar slab is achievable at all -- its negation is the truly-polar Tasker III verdict), whether a point flip survives adding vacuum (has_symmorphic_flip_op_by_bulk), and the full crystallographic surface-polarity class of Hinuma et al. (hinuma_polarity_type: polar / nonpolar_A / nonpolar_B / nonpolar_C, mapping onto the three Tasker ionic-surface types). pymatgen returns only a per-slab is_symmetric() bit and no orientation-level polarity classification (see Surface polarity (Hinuma type)).

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 + τ (g is a cleave's position along the normal, +/− keeps/flips the normal direction, and τ is the operation's translation along the normal). 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.

Terminology: termination vs. face

Throughout SlabTerminator, a termination means one symmetrically distinct cleave (a distinct way to cut the crystal along the Miller plane), and there is exactly one termination (one as-cut slab) per orbit of candidate cleavage gaps. get_unique_terminations(), n_unique_terminations, and termination_id all count and label cleaves in this sense.

A single cleave produces a slab with two faces (top and bottom). When the crystal has a flipping operation relating them the slab is nonpolar and both faces are equivalent; when it does not, the slab is polar and its two faces are inequivalent surface terminations in the surface-science sense. So one termination (one cleave) can expose two distinct faces, which is exactly why a polar slab yields two symmetrized children (one per face) and why include_flipped=True emits the opposite-face counterpart of each polar slab. To keep the vocabulary unambiguous, "termination" in this codebase always denotes the cleave, never a single face; the face notion appears only through include_flipped and the flipped field.

Benchmarking script details

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 for plain generation in scripts/benchmark_<date>.csv and for symmetrized generation in scripts/benchmark_symmetrize_<date>.csv, both date-stamped and co-measured in one interleaved run) for a correctness + speed comparison against the old UniqueSlabsGenerator and pymatgen's SlabGenerator (based on StructureMatcher).

By default the benchmark script recomputes all timings; the --reuse flag pulls whole method groups from the existing CSVs instead of re-measuring them, spanning both datasets: new/old/pmg (plain generation: SlabTerminator, old fingerprint, pymatgen) and new_sym/auto_sym/pmg_sym (symmetrized generation: ST at tol=0.1, ST at tol=auto, pymatgen). For example --reuse pmg,pmg_sym re-times the fast ST methods fresh while keeping both slow, cached pymatgen numbers. Each group keeps its own "measured at" timestamp in its meta sidecar. Using the flag --reuse all produces only the report. scripts/benchmark_analysis.py reads those CSVs (the most recent date on disk by default, or a specific run via --date YYYY-MM-DD) to plot the speedup-vs-size analysis (scripts/benchmark_speedup_vs_size_<date>.png, scripts/benchmark_symmetrize_speedup_vs_size_<date>.png, and the combined figure above) and writes scripts/benchmark_summary_<date>.md and scripts/benchmark_symmetrize_summary_<date>.md, each stamped with the analyzed run's date.

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

All surfaces of one material

To enumerate every symmetrically distinct surface of a bulk material, use the SlabTerminator.for_miller_family classmethod. It runs the bulk symmetry analysis once and hands the resulting operations to each per-facet SlabTerminator, so a whole family costs one spglib call instead of one per Miller index — the same reuse the batch pipeline does, in a lightweight form for library/interactive use:

from slabterminator.core import SlabTerminator

# One representative per distinct family up to |h|,|k|,|l| = 2 (e.g. (100), (110),
# (111), ... for cubic). Each yielded generator is ready to build slabs.
for miller, gen in SlabTerminator.for_miller_family(structure, max_miller_index=2):
    result = gen.get_unique_slabs(vacuum_size=15.0, min_slab_thickness=10.0)
    print(miller, result.properties.n_unique_terminations)

Construction only analyzes terminations; the vacuum/thickness/symmetrize knobs stay on the per-generator get_unique_slabs(...) call, so you keep full control (or skip a facet) per Miller index. Any keyword SlabTerminator accepts (tol, slab_symprec, output_format, ...) can be passed through — symprec is forwarded to both the shared bulk analysis and every generator, so pass it here rather than in the loop. For finished slab records across many materials in parallel, use slabterminator.pipeline / batch instead.

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, slab_symmetry_method="direct", output_format="pmg")

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, for a disordered (partially occupied) structure — the analysis needs a single species per site, so order the sites first — and for slab_symprec < symprec (see below).

symprec is the crystallographic tolerance on the clean bulk cell (it drives flip-op detection, the Hinuma class and orbit reduction); slab_symprec (defaulting to symprec) verifies the noisier with-vacuum slab face-symmetry during nonpolar trimming. They stay independent knobs, but slab_symprec must be the equal-or-coarser one — slab_symprec < symprec raises ValueError. This guarantees has_flip_op_by_bulk and the with-vacuum achievability verdict cannot diverge as a tolerance artifact (a flip op detected at a loose symprec but rejected by a tighter slab verification).

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.

slab_symmetry_method selects the backend for the with-vacuum face-symmetry check (is a built slab's two faces equivalent). The default "direct" applies the bulk flipping operations straight to the slab's atoms and tests self-coincidence in numpy — comparing only atom pairs within tolerance along the surface normal (roughly O(N × atoms-per-z-layer), at worst O(N²)) per flip operation, no per-slab spglib call. Because it works on the finite atom set and the surface {a, b} lattice only, its verdict is independent of both vacuum_size and max_normal_search; in particular it correctly detects a symmorphic (point) flip that SpacegroupAnalyzer can miss when the oriented cell's c vector is oblique. "spglib" runs a full SpacegroupAnalyzer on the built slab instead and is retained as a ground-truth oracle for benchmarking and testing (the two agree once spglib is given a non-oblique cell). Switching to "spglib" also brings back the low-vacuum_size warning, since only that backend can be fooled by periodic images across a thin gap.

output_format selects what each emitted SlabEntry.slab is: "pmg" (the default) a pymatgen Structure, or "dict" an as_dict()-shaped dict built straight from the cleave arrays — skipping both the per-slab Structure construction and the Structure.as_dict() serialization, which pays off when persisting large datasets (the SLURM build script sets it). A "dict" slab reconstructs to the same Structure via Structure.from_dict. "dict" requires the default slab_symmetry_method="direct" (the "spglib" backend needs a real Structure), so combining it with "spglib" raises ValueError.

• 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), symmetric_by_bulk (cheap pre-vacuum face-symmetry estimate), and termination_key (the cross-run-stable content key described under slabs above). 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.
  • include_flipped (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).
  • symmetrize (default False): for each polar slab, also emit up to two symmetrized slabs (usually nonstoichiometric) — one made nonpolar by trimming the top face, one by trimming the bottom — when the crystal has a flipping symmetry operation. The equivalent of pymatgen's symmetrize flag, but done in closed form: the slab is trimmed straight to a crystal flip center (each candidate verified with the with-vacuum symmetry check) instead of removing one atom at a time. Trimming removes the fewest atoms that still meet min_slab_thickness and keep the trimmed slab's fixed face equal to the parent's as-cut face (so a child reproduces the parent surface rather than exposing a deeper termination). Accepts a bool or a mode string controlling which polar slabs are emitted alongside the children: False/None off; True/"add" keeps every as-cut slab and flags a symmetrized parent has_symmetrized_children=True so it can be filtered out later; "drop_parents" drops a parent (and its flip) once it produced children but keeps a polar termination that could not be symmetrized (Tasker III); "drop_all_polar" keeps only nonpolar slabs — children plus any termination nonpolar as-cut — matching pymatgen's get_slabs(symmetrize=True). Under the drop modes a dropped slab's Structure is never built.

  The returned UniqueSlabsResult is a NamedTuple:

  • properties: n_unique_terminations, surface_area, two symmetry flags whose names advertise how they are determined — has_flip_op_by_bulk (a flipping projected operation exists for this crystal × Miller; this is exactly the "a nonpolar slab is achievable at all" verdict, so its negation False is the truly-polar, Tasker type III verdict), and has_symmorphic_flip_op_by_bulk (a point flip op exists — inversion / mirror ⊥ normal / in-plane 2-fold — an exact certificate that a nonpolar termination or trimming exists) — plus hinuma_polarity_type, the crystallographic surface-polarity class (see Surface polarity). They satisfy has_symmorphic_flip_op_by_bulk ⟹ has_flip_op_by_bulk.
  • 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, or an as_dict()-shaped dict when output_format="dict"), nsites and thickness_A (the slab's site count and occupied along-normal thickness in Angstrom, read straight from the arrays so they need no inspection of the slab payload), three classifier fields — flipped (bool, the turned-over counterpart), has_symmetrized_children (bool, a polar parent that produced symmetrized children) and symmetrized ('top', 'bottom' or None) — window_top (the top cleave for a symmetrized slab, else None) and boundary_crop ('drop'/'keep', the window's boundary treatment — pass it to regenerate_slabs(..., boundary_crops=[...]) to reproduce a parent-faithful kept-boundary child exactly), and off_stoichiometry (a {element: surplus} dict giving the atoms in excess of the largest whole bulk-formula multiple the slab contains — empty for every stoichiometric slab, including any symmetrized slab whose trim happened to preserve the bulk ratio; non-empty only for symmetrized ones), and termination_id (the 0-based index of the source termination — an as-cut parent, its flip, and its symmetrized children share one termination_id, so a symmetrized child is matched to its parent by equal termination_id, the parent being the entry with that id whose flipped is False and symmetrized is None). Two content-based join keys accompany the ordinal: termination_key — a family key (a 2-layer per-face composition + interlayer-spacing profile, e.g. O1@0.51|O1;Ti1@0.91|O1) shared by the parent/flip/children like termination_id but stable across runs/platforms at a fixed tol, so datasets join on it even when an unrelated termination is added or removed — and slab_key — a per-slab identity ({termination_key}#{role}:n{nsites}:{own_face}) that is instead distinct for the parent, flip and each child, so a specific sibling can be joined across runs. See the field docstrings for the two documented caveats (termination_key is not fully tol-invariant and is in-plane-degenerate for ~5% of high-symmetry pairs). Each SlabEntry also carries diagnostics, a (usually empty) tuple of per-slab provenance codes (see Diagnostics).

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

Rebuilds the final slabs 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. Pass flipped (the parallel SlabEntry.flipped values) to reproduce turned-over entries and windows (the window_top values) to reproduce symmetrized entries. Returns pymatgen Structures by default, or as_dict()-shaped dicts with output_format="dict".

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.

With symmetrize=True, tol="auto" picks the plateau of the full output — the terminations plus the symmetrized children each polar termination produces — rather than the termination count alone. Whether a termination is polar and how many symmetric trimmings it admits are properties of the cleave position (independent of tol, and of min_slab_thickness, vacuum, and the output cell), so the child-aware scan runs on the minimal analysis cell and memoizes each cleave's contribution across the swept tolerances. It is resolved lazily by get_unique_slabs(symmetrize=True) for that build alone and reported as settings.tol; gen.tol keeps the termination-count "auto" value, so the same object can report a different (child-aware) tolerance for a symmetrize=True build than for a symmetrize=False one.

Diagnostics

A build records any provenance concerns it hits as a closed set of string codes, surfaced two ways at once: a durable field — SlabProperties.diagnostics (orientation-level) and SlabEntry.diagnostics (per-slab), unioned onto SlabRecord.diagnostics and the SLURM diagnostics column — and a live warnings.warn at the origin. The field is empty on a clean build, so a finished dataset can be filtered for "which rows had a concern" rather than scraping logs. There are currently two codes:

• tol_ambiguous — the effective tol="auto" scan found no stable plateau (the ambiguous case above), so the unique count is genuinely resolution-dependent and 0.1 was kept. Orientation-level (on SlabProperties).
• trim_crop_unresolved — a symmetrized child's fixed face reproduced neither the drop nor the keep crop of the parent surface layer in the output cell, so the analysis-cell crop was kept unverified. A believed-invariant defensive case (not seen on the index ≤ 3 dataset); per-slab (on SlabEntry).

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)

SlabTerminator.for_miller_family (see All surfaces of one material) wraps exactly this loop when you just want every distinct facet of a material.

Symmetry flags

Each boolean symmetry tag is named for how it is determined. Orientation-level flags belong to the (crystal, Miller index) pair (on SlabProperties, mirrored onto every SlabRecord); termination-level flags describe one termination (cleave) or its slab.

The From column has exactly two primitives, plus two combinations of them:

• bulk ops — algebra on the bulk crystal's projected symmetry operations; no slab is built.
• verified — decided by building the with-vacuum slab and checking it (default direct backend, or spglib — see slab_symmetry_method).
• bulk ops, else verified — settled by algebra where it can be, building a slab to verify only the ambiguous case.
• bulk ops + stoichiometry — algebra plus a per-layer composition check: whether each atomic (hkl) layer, on its own, already has the bulk element ratio.

Flag	Level	From	Meaning
has_flip_op_by_bulk	orientation	bulk ops	a flipping op exists (normal → −normal). This is exactly the "a nonpolar slab is obtainable at all" verdict; False is truly polar (Tasker III, unreconstructable): no nonpolar termination is achievable (not merely "has a polar termination").
has_symmorphic_flip_op_by_bulk	orientation	bulk ops	a point flip op exists (intrinsic glide/screw translation ≈ 0) — exact certificate that a nonpolar termination or trimming exists.
hinuma_polarity_type	orientation	bulk ops + stoichiometry	crystallographic surface-polarity class: polar / nonpolar_A / nonpolar_B / nonpolar_C (see Surface polarity).
symmetric_by_bulk	termination	bulk ops	a bulk flip fixes this cleave; over-estimates face symmetry.
is_symmetric_with_vacuum	termination	verified	the built slab's two faces are genuinely equivalent (nonpolar as-cut).
has_symmetrized_children	termination	verified	this polar parent produced verified symmetrized children (symmetrize=True only).

Relationships

The names are chosen so these implications read directly — the left-hand flag implies the right-hand one for that crystal × Miller (and where the left is a per-termination flag, it implies its orientation-level counterpart):

is_symmetric_with_vacuum        ⟹  symmetric_by_bulk  ⟹  has_flip_op_by_bulk
has_symmorphic_flip_op_by_bulk  ⟹  has_flip_op_by_bulk
has_symmetrized_children        ⟹  has_flip_op_by_bulk

has_flip_op_by_bulk is itself the "a nonpolar slab is achievable at all" verdict: a bulk flip always admits a nonpolar trimming (the trim re-centers the slab on the flip center, so even a non-symmorphic glide/screw flip survives adding vacuum), and no flip op means no termination can ever be nonpolar. __init__ enforces slab_symprec >= symprec (bulk detection at the finer symprec, with-vacuum verification at the equal-or-coarser slab_symprec), so the two verdicts cannot diverge as a tolerance artifact.

Contrapositively, negations cascade the other way; existential antecedents ("some termination") become universal ("every termination", marked ∀ below):

# No flipping op at all — the strongest polar condition (truly polar, no nonpolar slab
# obtainable), so nothing is symmetric:
¬has_flip_op_by_bulk  ⟹  ¬has_symmorphic_flip_op_by_bulk
                      ⟹  ¬symmetric_by_bulk         (∀ terminations)
                      ⟹  ¬is_symmetric_with_vacuum  (∀ terminations — none nonpolar as-cut)
                      ⟹  ¬has_symmetrized_children  (∀ entries — nothing symmetrized)

# Per-termination short-circuit (skips the with-vacuum check in the code):
¬symmetric_by_bulk    ⟹  ¬is_symmetric_with_vacuum  (same termination)

The _by_bulk flags are cheap over-estimates, and the arrows do not reverse — the gaps are themselves meaningful:

• has_flip_op_by_bulk ∧ ¬has_symmorphic_flip_op_by_bulk — flip ops exist but are all non-symmorphic (glide/screw), so none is a point op that is nonpolar as-cut. This does not make the orientation polar: trimming re-centers the slab on the flip center, so a nonpolar slab is still obtainable (has_flip_op_by_bulk stays True, and a nonpolar slab is achievable, confirmed by a with-vacuum trimming check).
• symmetric_by_bulk ∧ ¬is_symmetric_with_vacuum — a termination whose bulk flip is broken by the vacuum.
• has_flip_op_by_bulk ⇏ has_symmetrized_children — an orientation is nonpolar via an as-cut nonpolar termination, or a certified trimming too thin to emit. Hence has_flip_op_by_bulk is the exact, build-parameter-independent achievability verdict; the emitted is_symmetric_with_vacuum / symmetrized slabs are a thickness-dependent view that can under-report.

Surface polarity (Hinuma type)

hinuma_polarity_type labels the (crystal, Miller index) with the crystallographic surface-polarity class of Hinuma et al. (Comput. Mater. Sci. 113 (2016) 221), computed purely from symmetry and stoichiometry (no charge, oxidation state, or dipole is used — matching the paper). It is one of four values:

Type	Condition	Example	~ Tasker
polar	no flip op exists (¬has_flip_op_by_bulk)	zincblende (111), wurtzite (0001)	III (polar instability)
nonpolar_A	flip op exists and every atomic (hkl) plane is on its own stoichiometric	zincblende (110), rocksalt (100)	I
nonpolar_B	flip op exists, not type A, but at least one simply-cleaved slab is both nonpolar and stoichiometric	fluorite (111)	II
nonpolar_C	flip op exists but no simply-cleaved slab is both nonpolar and stoichiometric — a nonpolar stoichiometric slab needs reconstruction	fluorite (100), rocksalt (111)	III (reconstruction)

Two invariants tie it to the flags above:

hinuma_polarity_type == "polar"                       ⟺  ¬has_flip_op_by_bulk
hinuma_polarity_type ∈ {"nonpolar_A", "nonpolar_B", "nonpolar_C"}  ⟹  has_flip_op_by_bulk

nonpolar_C vs polar — both map to "Tasker III", but they are distinct cases. The difference is where the crystal's flip symmetry lives:

• nonpolar_C: a flip symmetry does exist, but it sits on atomic planes rather than in the gaps between them. So a simply-cleaved stoichiometric slab is polar only because the cut is off-center; a stoichiometry-preserving reconstruction (e.g. the octopolar removal of half the atoms on both boundary planes) restores a nonpolar stoichiometric slab. Example: rocksalt (111) — alternating pure-cation / pure-anion planes with inversion centers on the atoms; centring a slab on a plane makes the faces match but adds an extra plane (non-stoichiometric), and cleaving in a gap keeps stoichiometry but leaves unlike faces — you can have one or the other, not both, until you reconstruct.
• polar: no flip symmetry exists anywhere. The two faces can never be made equivalent by any cut or stoichiometric reconstruction; only electronic/chemical compensation (charge redistribution, adsorbates, defects) resolves the dipole. Example: wurtzite (0001) — c is a polar axis, so there is no mirror or inversion ⊥ normal to exploit.

The discriminator is exactly has_flip_op_by_bulk: False → polar; True with no stoichiometric nonpolar simple cleave → nonpolar_C.

The A/B/C ≈ Tasker 1/2/3 correspondence holds barring exceptions: Tasker's scheme is defined on per-layer formal charge, which is deliberately not computed here, so charge-neutral-but-not-stoichiometric layers differ — e.g. SrTiO₃ (001) is Hinuma nonpolar_C (SrO / TiO₂ planes aren't individually stoichiometric) but Tasker type 1 (those planes are formally charge-neutral). A future formal-charge/oxidation-state per-layer analysis (TODO in the code) would enable emitting true Tasker types; until then no Tasker column is written, to avoid implying a rigor the geometry alone can't give.

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.slab_id, rec.symmetrized or ("flipped" if rec.flipped else "as_cut"),
          "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 (an as_dict()-shaped dict under
    # SlabGenConfig(output_format="dict")); regenerate_slabs([rec.shift], 0.0,
    # rec.oriented_unit_cell, flipped=[rec.flipped], windows=[rec.window_top]) rebuilds
    # the no-vacuum slab.

Each SlabRecord carries the built slab, its shift, oriented_unit_cell, and the per-slab and aggregate properties (is_symmetric_with_vacuum, has_flip_op_by_bulk, has_symmorphic_flip_op_by_bulk, hinuma_polarity_type, surface_area, n_unique_terminations, the resolved tol and max_normal_search, flipped, has_symmetrized_children, symmetrized, window_top, off_stoichiometry, termination_id (the source termination's index, linking children to their parent), …). The no-vacuum slab is not stored: the oriented cell plus the shift (and flipped / window_top) reconstruct it exactly via regenerate_slabs. Set SlabGenConfig(symmetrize=True) to additionally emit the nonstoichiometric symmetrized slabs described above.

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 23,400 tests, most of them parametrized across the fixtures below (39 fixture structures in tests/test-cifs/; the directory also holds two un-parametrized CIFs loaded directly by individual regression tests: TiNi (mp-1179013), pulled for a count-mismatch investigation, and PaAs (mp-11106), a symmetrized-child as-built-count-guard case).

The core set is 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: 31 Materials Project conventional standard cells plus one synthetic polar structure (synthetic_polar_Pna21.cif) for the point group mm2. Every fixture in this set is run through the same parametrized sweeps. Two extra named fixtures are included beyond the point-group survey as regression cases: Ce2NiGe3 (mp-1102475), a glide / rich-symmetry cell whose reversal-canonical lower bound over-counts, and Tb5Ir3 (mp-1188649), a P-hexagonal metal whose (1,3,3) surface exercises the symmetrized-window max_normal_search invariance.

The higher-index sweep (every symmetrically distinct Miller index up to 3, for every fixture) is validated without hand-curated counts: each case is checked against independent, orbit-free bounds — a hard gap-count ceiling, a symmetry-gated lower floor, and shift-free count brackets — backed by an independently re-derived projected-orbit oracle for the screw / sub-period regime those bounds cannot pin.

Five further MP reference cells back the surface-polarity (hinuma_polarity_type) tests against the paper's own published classifications — MgO rocksalt ((100) type A, (111) type C), SrTiO₃ perovskite ((001) type C, the Hinuma-vs-Tasker exception), and the paper's three worked type-B examples BeSO₄, high-pressure AgI, and FeSe₂; wurtzite AlN ((0001) polar) reuses the point-group cell. These are kept out of the point-group survey so it stays one structure per point group.

References

The hinuma_polarity_type surface-polarity classification (polar / nonpolar_A / nonpolar_B / nonpolar_C) implements the crystallographic scheme of:

Y. Hinuma, Y. Kumagai, F. Oba, I. Tanaka, "Categorization of surface polarity from a crystallographic approach", Computational Materials Science 113 (2016) 221–230. doi:10.1016/j.commatsci.2015.11.042

The A/B/C classes map (barring formal-charge exceptions) onto the three ionic-surface types of:

P. W. Tasker, "The stability of ionic crystal surfaces", Journal of Physics C: Solid State Physics 12 (1979) 4977–4984. doi:10.1088/0022-3719/12/22/036