Metadata-Version: 2.5
Name: slabterminator
Version: 0.5.2
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) cases across 1000 materials (run with
  [`scripts/benchmark.py`](scripts/benchmark.py)), `SlabTerminator` is, per case, roughly
  **15× (index 1) to 89× (index 3) faster** than `SlabGenerator` when building slabs — a
  **geometric-mean ~60× across the full dataset** (median ~72×; the tail-weighted
  whole-workload aggregate is ~95×) — and **~109× 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.80×) 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`](scripts/benchmark.py) run), the **per-case**
  generation speedup over pymatgen is even more pronounced: geometric-mean roughly **17×
  (index 1) to 160× (index 3)** per case (**~94× across the full dataset**; the tail-weighted
  aggregate is ~217×). 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 advantage 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**).

  <div align="center"><img src="scripts/benchmark_combined_speedup_vs_size_2026-08-17.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. 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 aggregate ~217× 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`](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.meta.md`](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](#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)](#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`](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
```

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

### `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)
```

## 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](#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 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:

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

```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 15,900 tests, most of them parametrized across the
fixtures below (38 fixture structures in [`tests/test-cifs/`](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](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)
