Metadata-Version: 2.4
Name: slabterminator
Version: 0.3.1
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">
  <img alt="SlabTerminator Logo" src="logo.svg" width="360"><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.** On our benchmark of 1352 (structure, Miller) cases,
  `SlabTerminator` is roughly **40 to 131× faster** than `SlabGenerator` when building
  slabs (**111× on average**) and **~320× faster on average** when only counting terminations,
  because it skips both the per-cleave slab construction and the pairwise `StructureMatcher`
  deduplication. In `batch` mode there is an additional **~26%** on top as the bulk symmetry is
  computed once per material and reused across all its Miller indices.

  <div align="center"><img src="scripts/benchmark_speedup_vs_size.png" width="70%" alt="Speedup of SlabTerminator over pymatgen SlabGenerator as a function of problem size and Miller index"></div>

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

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

## How it works

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

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

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

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

## A note on AI-assisted development

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

## Installation

Requires Python ≥ 3.12. The package is on [PyPI](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)`

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

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

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

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

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

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

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

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

### Choosing the tolerance automatically

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

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

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

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

## 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.term_index, rec.face,
          "polar" if not rec.is_symmetric_with_vacuum else "nonpolar",
          rec.top_layer_composition, "/", rec.bottom_layer_composition)
    # rec.slab is the with-vacuum pymatgen Structure; rec.oriented_unit_cell + rec.shift
    # reconstruct the no-vacuum slab via regenerate_slabs(vacuum_size=0.0).
```

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

### `slabterminator.batch` — many materials in parallel

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

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