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) cases across 1000 materials (run with scripts/benchmark.py), SlabTerminator is on average roughly 27× (index 1) to 128× (index 3) faster per case than SlabGenerator when building slabs (~101× across the full dataset) and ~216× faster when only counting terminations, because it skips both the per-cleave slab construction and the pairwise StructureMatcher deduplication. In batch mode there is a further ~1.7× on top (up to 2.85×) as the bulk symmetry is computed once per material and reused across all its Miller indices.

  When the slabs are being symmetrized (symmetrize=True, co-measured in the same scripts/benchmark.py run), the per-case generation speedup over pymatgen is even more pronounced: on average roughly 31× (index 1) to 272× (index 3) per case (~198× across the full dataset). 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 centre in closed form. That ratio is dominated by pymatgen's tail; end-to-end on a full dataset build the realistic advantage is smaller, about 33× per slab / 20× per material, because SlabTerminator also emits ~1.7× more (verified-nonpolar) surfaces (see the note below). pymatgen additionally trims each polar slab atom-by-atom and re-runs its symmetry analysis (detailed speedup breakdown, see right plot).

  

  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 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. An independent HPC run generating slabs up to a Miller index of 3 for 1000 materials in full batch mode (10 workers on Northeastern University's Explorer HPC Cluster, build + serialization + JSONL I/O) took 112.0 s with SlabTerminator versus 10121.4 s with pymatgen — a ~90× speedup on real production hardware at a matched fixed layer tolerance (tol = ftol = 0.1). With symmetrize=True the same build took 580.3 s versus 11446.1 s (~20×). That end-to-end symmetrize ratio is lower than the per-case ~198× above mainly because SlabTerminator emits (and serializes) ~1.7× more, more-complete slabs — the polar parents and unsymmetrizable Tasker-III terminations pymatgen drops — and the batch pays a per-slab serialization/write tail (shared by both generators) that the generation-only per-case figure excludes. 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.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), whether a point flip survives adding vacuum (has_symmorphic_flip_op_by_bulk), whether any nonpolar slab is achievable at all (nonpolar_slab_achievable, the truly-polar Tasker III verdict), 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 + τ. 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

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

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

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 — 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), 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.
  • 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 centre (each candidate verified with the with-vacuum symmetry check) instead of removing one atom at a time. Trimming removes the fewest atoms while still meeting min_slab_thickness. 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, three symmetry flags whose names advertise how they are determined — has_flip_op_by_bulk (a flipping projected operation exists for this crystal × Miller; necessary but not sufficient for a nonpolar surface), 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), and nonpolar_slab_achievable (whether a nonpolar slab is obtainable at all; False is the truly-polar, Tasker type III verdict, computed on every run regardless of symmetrize) — plus hinuma_polarity_type, the crystallographic surface-polarity class (see Surface polarity). They satisfy has_symmorphic_flip_op_by_bulk ⟹ nonpolar_slab_achievable ⟹ 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 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).

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.

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)

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 or its slab. _by_bulk flags are algebra on the bulk crystal's projected symmetry ops (no slab built); the rest are verified from the built, with-vacuum slab (the default direct backend, or spglib — see slab_symmetry_method) or record its outcome.

Flag	Level	From	Meaning
has_flip_op_by_bulk	orientation	bulk ops	a flipping op exists (normal → −normal). Necessary, not sufficient, for a nonpolar surface.
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.
nonpolar_slab_achievable	orientation	bulk ops + verified	a nonpolar slab is obtainable at all. False is truly polar (Tasker III, unreconstructable): no nonpolar termination is achievable (not merely "has a polar termination").
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 outcome	this polar parent produced verified symmetrized children (symmetrize=True only).

Relationships

The names are chosen so these implications read directly (a per-termination flag on the left implies the orientation-level flag on the right for that crystal × Miller):

is_symmetric_with_vacuum        ⟹  symmetric_by_bulk  ⟹  has_flip_op_by_bulk
has_symmorphic_flip_op_by_bulk  ⟹  nonpolar_slab_achievable  ⟹  has_flip_op_by_bulk
has_symmetrized_children        ⟹  nonpolar_slab_achievable
is_symmetric_with_vacuum        ⟹  nonpolar_slab_achievable

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:
¬has_flip_op_by_bulk  ⟹  ¬nonpolar_slab_achievable                              (definitional)
                      ⟹  ¬has_symmorphic_flip_op_by_bulk
                      ⟹  ¬symmetric_by_bulk         (∀ terminations)
                      ⟹  ¬is_symmetric_with_vacuum  (∀ terminations)
                      ⟹  ¬has_symmetrized_children  (∀ entries)

# Truly polar (no nonpolar slab obtainable) ⇒ nothing is symmetric:
¬nonpolar_slab_achievable  ⟹  ¬has_symmorphic_flip_op_by_bulk
                           ⟹  ¬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 the vacuum may break them; nonpolar_slab_achievable resolves this case with a with-vacuum trimming check.
• symmetric_by_bulk ∧ ¬is_symmetric_with_vacuum — a termination whose bulk flip is broken by the vacuum.
• ¬nonpolar_slab_achievable ⇏ ¬has_flip_op_by_bulk — a truly polar orientation can still have flip ops, if all are non-symmorphic and all break under vacuum.
• nonpolar_slab_achievable ⇏ has_symmetrized_children — an orientation is nonpolar via an as-cut nonpolar termination, or a certified trimming too thin to emit. Hence nonpolar_slab_achievable is the exact, build-parameter-independent 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 a 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_slab_achievable

nonpolar_C vs polar — both are "Tasker III", but they are not the same problem. 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-centre; 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 centres 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.

hinuma_polarity_type (crystallographic class) and nonpolar_slab_achievable (physical: can this tool actually build a nonpolar slab, vacuum-aware) answer different questions and agree except in one corner: when the only flip ops are non-symmorphic (glide/screw), the infinite crystal has a flip op (so Hinuma reports nonpolar_C, reconstructable in principle), yet the vacuum may break every such op so that no nonpolar slab is actually buildable (nonpolar_slab_achievable == False). There, the two disagree by design.

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, nonpolar_slab_achievable, 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 15,900 tests, most of them parametrized across the fixtures below (38 fixture structures in tests/test-cifs/; the directory also holds a handful of un-parametrized CIFs pulled for count-mismatch investigations).

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. One extra named fixture, Ce2NiGe3 (mp-1102475), is included beyond the point-group survey as a regression case for the reversal-canonical bracket — a glide / rich-symmetry cell whose lower fingerprint bound over-counts.

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