Metadata-Version: 2.5
Name: slabterminator
Version: 0.6.0
Summary: Enumerate the symmetrically unique slab terminations of a bulk crystal for a given Miller index from crystal symmetry, and generate slabs.
Project-URL: Homepage, https://github.com/d2r2group/slabterminator
Project-URL: Repository, https://github.com/d2r2group/slabterminator
Author-email: Peter Schindler <p.schindler@northeastern.edu>
License-Expression: MIT
License-File: LICENSE.md
Keywords: Miller index,crystallography,materials science,pymatgen,slab,surface science,surface termination,symmetry
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.12
Requires-Dist: pymatgen>=2026.5.4
Description-Content-Type: text/markdown

<div align="center">
  <picture>
    <source media="(prefers-color-scheme: dark)" srcset="logo.svg">
    <img alt="SlabTerminator Logo" src="logo-light.svg" width="360">
  </picture><br>
</div>

# SlabTerminator

![Python - Version](https://img.shields.io/pypi/pyversions/slabterminator)
[![PyPI - Version](https://img.shields.io/pypi/v/slabterminator?color=blue)](https://pypi.org/project/slabterminator)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

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](#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`](scripts/benchmark.py) and the full 
  list of mpids used for benchmarking is stored in `benchmark_HPC_<date>.meta.md`.

  <div align="center"><img src="scripts/benchmark_combined_speedup_vs_size_2026-08-18.png" width="100%" alt="Speedup of SlabTerminator over pymatgen SlabGenerator for symmetrize=False (left) and symmetrize=True (right), as a function of problem size and Miller index, with y-axes matched across the two"></div>

  *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`](scripts/build_slab_dataset_slurm.py)
  and [`scripts/build_slab_dataset_slurm_pymatgen.py`](scripts/build_slab_dataset_slurm_pymatgen.py),
  matched on everything but the generator (run-wide details in
  [`scripts/benchmark_HPC_<date>.meta.md`](scripts/benchmark_HPC_2026-08-18.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](#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)](#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`](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`](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](https://pypi.org/project/slabterminator):

```bash
# 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`](https://docs.astral.sh/uv/)
to install it with its dev dependencies:

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

## Quick start

```python
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](#slabterminatorpipeline-one-material-all-miller-indices) does, in a
lightweight form for library/interactive use:

```python
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`](#batch-generation-over-many-materials) 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](#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](#surface-polarity-hinuma-type)). 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](#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 `Structure`s 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):

```python
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):

```python
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:

```python
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](#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](#surface-polarity-hinuma-type)). |
| `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:

```python
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.

```python
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

```bash
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/`](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](https://doi.org/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](https://doi.org/10.1088/0022-3719/12/22/036)
