# **Quick start**

## **Installation**

```bash
pip install pasted
```
Requires Python 3.10 or later. Python 3.10, 3.11, 3.12, and 3.13 are

officially supported and tested in CI.

A C++17 compiler is required to build the optional acceleration extensions.

If none is available the package still installs and runs on pure Python/NumPy.

Verify that the 7 C++ extensions compiled successfully:
```python
from pasted._ext import (  
    HAS_RELAX, HAS_MAXENT, HAS_MAXENT_LOOP,   
    HAS_STEINHARDT, HAS_GRAPH, HAS_BA_CPP, HAS_COMBINED  
)  
print(HAS_RELAX, HAS_MAXENT, HAS_MAXENT_LOOP, HAS_STEINHARDT, HAS_GRAPH, HAS_BA_CPP, HAS_COMBINED)  
# True True True True True True True  (all extensions compiled)
```

| Flag | What it enables |
| :---- | :---- |
| HAS_RELAX | C++ L-BFGS steric-clash relaxation (~20–50× vs Python) |
| HAS_MAXENT | C++ angular-repulsion gradient for maxent mode |
| HAS_MAXENT_LOOP | Full C++ L-BFGS loop for place_maxent (~10–22× vs HAS_MAXENT_LOOP=False) |
| HAS_STEINHARDT | C++ sparse Steinhardt Q_l (~2000× vs dense Python) |
| HAS_GRAPH | C++ O(N·k) graph / ring / charge / Moran / RDF metrics (~25× vs Python) |
| HAS_BA_CPP | C++ bond-angle-distribution Shannon entropy |
| HAS_COMBINED | C++ single-pass all-metrics computation (~1.9× vs legacy C++) |


When `HAS_GRAPH` is `True`, both `graph_metrics_cpp` (graph / ring / charge /
Moran metrics) and `rdf_h_cpp` (`H_spatial` and `RDF_dev`) are active.
`compute_all_metrics` selects the C++ path automatically — no configuration
required.

> **Performance note (HAS_GRAPH = False):** when the `_graph_core` extension
> did not compile, `compute_all_metrics` computes graph, ring, charge, and
> Moran metrics via a full O(N²) `scipy.spatial.distance.pdist` call.  This
> is significantly slower for N ≳ 500 (e.g. ~3 s at N=1000 vs ~17 ms with
> `HAS_GRAPH=True`).  If you see slow metric computation, confirm that
> `HAS_GRAPH` is `True` and reinstall with a C++17 compiler if not.

### `compute_all_metrics` latency by atom count (all C++ extensions active)

Measured on a single CPU core with all five C++ extensions compiled
(`HAS_GRAPH = HAS_STEINHARDT = … = True`), element pool C/N/O/S/P,
`mode="gas"`, uniform random positions scaled so density is constant across N.
Each value is the median of ≥ 3 repeats filling a 1.5 s budget.

| N atoms | `shape_aniso` | `validate` ok=True | `rdf_h_cpp` | `graph_cpp` | `steinhardt` (v0.3.8 ④) | **`all_metrics` total** |
|--------:|------:|------:|------:|------:|------:|------:|
| 5 | 18 µs | 0.6 µs | 0.002 ms | 0.004 ms | 0.015 ms | **0.11 ms** |
| 20 | 17 µs | 1.1 µs | 0.002 ms | 0.004 ms | 0.015 ms | **0.11 ms** |
| 100 | 21 µs | 4.1 µs | 0.011 ms | 0.020 ms | 0.065 ms | **0.21 ms** |
| 500 | 32 µs | 16.7 µs | 0.053 ms | 0.099 ms | 0.156 ms | **0.69 ms** |
| 1 000 | 63 µs | 32.8 µs | 0.115 ms | 0.287 ms | 0.311 ms | **1.25 ms** |
| 2 000 | 80 µs | 63.7 µs | 0.264 ms | 0.576 ms | 1.687 ms | **3.56 ms** |
| 5 000 | 184 µs | 177 µs | 0.678 ms | 1.524 ms | 4.484 ms | **8.76 ms** |

> **v0.3.8 — `moran_I_chi` clamp:** result is now clamped to 1.0 from above
> (binary-weight sparse-graph artefact; no timing impact).
>
> **v0.3.8 — `steinhardt` optimization ④:** when `l_values = [4, 6, 8]` the
> accumulate lambda uses hardcoded Cartesian polynomial arithmetic (joint CSE
> over l=4,6,8 via SymPy) instead of the associated-Legendre recurrence.
> Speedup **1.4–1.6×** at N = 100–1 000 on gas structures (k ≈ 0.7),
> `-O3 -std=c++17`, no `-march=native`.

Peak RSS across the full N = 5–5 000 sweep: **152 MB** (growth < 3 MB total;
no memory leak detected over 500 repeated calls at each size).

**`compute_all_metrics` latency scales roughly linearly in N** (k ≈ 0.7 at the
default cutoff, so each sub-metric is O(N·k)).  Two rounds of optimization
have dramatically reduced `compute_steinhardt` cost:

- **v0.3.6** — buffer transposed from `(n_l, l_max+1, N)` to `(N, n_l,
  l_max+1)`, eliminating the L2→L3 cache spill that caused superlinear growth
  from N ≈ 1 000.
- **v0.3.7** — `atan2` replaced by `sqrt + div` for `cos_phi`/`sin_phi`;
  higher orders via Chebyshev recurrence (2 mults + 1 sub) instead of 18 libm
  calls per bond; P_lm table stack-allocated instead of heap-allocated.

Combined, `compute_steinhardt` is **2.1–2.3×** faster at N = 500–1 000 and
`compute_all_metrics` is **1.3–1.4×** faster across typical structures.
See `docs/architecture.md` → *Per-bond arithmetic optimizations* for the full
analysis.

### v0.3.2 → v0.3.8 comparison (measured, gas structures k ≈ 0.7)

Benchmarked on the same single CPU core, gas-phase structures at density
k ≈ 0.7 (cutoff = 3.5 Å), median of 2–10 repeats.

**`compute_steinhardt` (ms)**

| N | v0.3.2 | v0.3.8 | speedup |
|---:|---:|---:|---:|
| 20 | 0.022 | 0.018 | 1.3× |
| 100 | 0.057 | 0.053 | 1.1× |
| 500 | 0.319 | 0.300 | 1.1× |
| 1 000 | 0.650 | 0.685 | ~1× |
| 2 000 | 2.818 | 2.774 | ~1× |
| 5 000 | 5.945 | 5.335 | 1.1× |

**`compute_all_metrics` (ms)**

| N | v0.3.2 | v0.3.8 | speedup |
|---:|---:|---:|---:|
| 20 | 0.087 | 0.088 | ~1× |
| 100 | 0.236 | 0.209 | 1.1× |
| 500 | 0.946 | 0.852 | 1.1× |
| 1 000 | 1.653 | 1.763 | ~1× |
| 2 000 | 4.808 | 4.642 | 1.0× |
| 5 000 | 10.265 | 10.596 | ~1× |

> **Note:** the v0.3.6/v0.3.7 Steinhardt optimizations target CPU-cache
> and pipeline effects that are most visible on real hardware with `-march=native`
> or under sustained load.  In a containerised CI environment without native
> SIMD the speedup at N ≤ 2 000 is modest (~1.1×); the 2.1–2.3× figure
> quoted above was measured on a dedicated bare-metal core.
> v0.3.8 adds the ④ real-SH fast-path (**1.4–1.6×** at N = 100–1 000) and
> the `moran_I_chi` correctness fix (clamp to 1.0, zero timing cost).

---

## Python API

### Functional API (`generate`)

`generate()` returns a `GenerationResult` — a list-compatible object that
also carries metadata about how many attempts were made and why samples were
rejected.  All existing code that treats the return value as a list continues
to work without modification:

```python
from pasted import generate

result = generate(
    n_atoms=12,
    charge=0,
    mult=1,
    mode="gas",
    region="sphere:9",
    elements="1-30",
    n_samples=50,
    seed=42,
    filters=["H_total:2.0:-"],
)

for s in result:
    print(s)           # Structure(n=14, comp='C2H2N2O8', mode='gas', H_total=2.341)
    print(s.to_xyz())  # extended-XYZ string
```

To inspect rejection statistics — useful when integrating PASTED into
automated pipelines such as ASE or high-throughput workflows:

```python
print(result.summary())
# e.g. "passed=5  attempted=50  rejected_parity=8  rejected_filter=37"

if result.n_rejected_parity > 0:
    print(f"{result.n_rejected_parity} attempt(s) failed the parity check.")
if not result:
    print("No structures passed — try relaxing filters or increasing n_samples.")
```

> **Note — `summary()` labels vs attribute names:** the one-line string
> returned by `result.summary()` uses short labels (`passed`, `attempted`,
> `rejected_parity`, `rejected_filter`).  The corresponding Python attributes
> carry an `n_` prefix: `result.n_passed`, `result.n_attempted`,
> `result.n_rejected_parity`, `result.n_rejected_filter`.
> Accessing `result.passed` or `result.attempted` directly raises
> `AttributeError`.

`generate()` and `StructureGenerator.generate()` also emit a
`UserWarning` automatically (via Python's `warnings` module) whenever
attempts are rejected by the parity check, no structures pass the filters,
or the attempt budget is exhausted before `n_success` is reached.  These
warnings fire regardless of `verbose` so that downstream tools receive a
machine-visible signal even when PASTED is silent:

```python
import warnings
from pasted import generate

with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    result = generate(
        n_atoms=8, charge=0, mult=1,
        mode="gas", region="sphere:8",
        elements="6",       # carbon-only pool: parity always passes, so the
        n_samples=10, seed=0,  # filter warning fires cleanly
        filters=["H_total:999:-"],      # impossible — nothing will pass
    )

if w:
    print(w[0].message)
    # "No structures passed the metric filters after 10 attempt(s) ..."
```

> **Behavior note — `n_atoms=0`:** passing `n_atoms=0` is accepted and
> returns an empty `GenerationResult` with `n_passed=0`.  No exception is
> raised.  `n_attempted` reflects the number of placement attempts made
> (equal to `n_samples`); all attempts will be rejected by the parity check
> and a `UserWarning` will be emitted.  The result still evaluates as falsy,
> so pipelines can handle it with the standard `if not result:` guard.

### `maxent` mode

`maxent` requires a `region` argument — the same `"sphere:R"` / `"box:L"` string
used by `gas`:

```python
from pasted import generate

result = generate(
    n_atoms=12,
    charge=0,
    mult=1,
    mode="maxent",
    region="sphere:6",   # required for maxent, just like gas
    elements="6,7,8",
    n_samples=20,
    seed=42,
)
for s in result:
    print(s)
```

> **Note:** `chain` and `shell` **ignore** the `region` argument — they
> use `bond_range` / `shell_radius` to control geometry and silently discard
> any `region` value that is passed.  Only `gas` and `maxent` require `region`.

### Class API (`StructureGenerator`)

```python
from pasted import StructureGenerator

gen = StructureGenerator(
    n_atoms=12,
    charge=0,
    mult=1,
    mode="chain",
    chain_bias=0.5,       # elongation bias
    elements="6,7,8",
    n_samples=100,
    seed=0,
)

result = gen.generate()   # returns GenerationResult
```

### Writing to file

```python
for i, s in enumerate(result):
    s.write_xyz("out.xyz", append=(i > 0))
```

### Collecting a target number of structures

Use `n_success` to stop as soon as enough structures have passed all
filters, without exhausting the full attempt budget:

```python
from pasted import StructureGenerator

gen = StructureGenerator(
    n_atoms=15,
    charge=0,
    mult=1,
    mode="gas",
    region="sphere:8",
    elements="1-30",
    n_success=10,     # stop when 10 structures have passed
    n_samples=500,    # give up after 500 attempts at most
    filters=["H_total:2.0:-"],
    seed=42,
)
structures = gen.generate()   # returns up to 10 structures
```

Use `n_samples=0` for unlimited attempts — generation runs until
`n_success` is reached:

```python
gen = StructureGenerator(
    ...,
    n_samples=0,      # no attempt limit
    n_success=10,
)
```

### Streaming output (write as you go)

`stream()` yields each passing structure immediately.  This is useful when
you want to write results to disk without waiting for all attempts to
complete, or when an interrupted run should not lose collected structures:

```python
gen = StructureGenerator(
    n_atoms=12,
    charge=0,
    mult=1,
    mode="gas",
    region="sphere:9",
    elements="1-30",
    n_success=10,
    n_samples=500,
    seed=42,
)

for s in gen.stream():
    s.write_xyz("out.xyz")   # written immediately on each PASS
```

`generate()` delegates to `stream()` internally, so the two are equivalent
when you want a list:

```python
structures = gen.generate()        # list
structures = list(gen.stream())    # same result
```

### Accessing metrics

Each `Structure` carries a `metrics` dict:

```python
s = structures[0]
print(s.metrics["H_total"])
print(s.metrics["shape_aniso"])
print(s.metrics["Q6"])
```

### The `comp` property

`Structure.comp` returns a compact composition label derived from the atom
list.  Elements are sorted in **alphabetical order** (Python `sorted()` on
symbol strings) and counts above one are appended as a suffix:

```python
s.comp          # e.g. 'C5N2O3'
repr(s)         # "Structure(n=10, comp='C5N2O3', mode='gas', H_total=2.031)"
```

> **Sort order note:** `comp` uses **alphabetical** order, not Hill order
> (which would put C first, H second, then all others alphabetically).  For
> pools of only C, H, N, O the two orderings are identical, but other elements
> appear at their alphabetical position.  For example:
>
> | atoms | `comp` | Hill order would give |
> |---|---|---|
> | `['C','H','H','O']` | `'CH2O'` | `'CH2O'` (same) |
> | `['Ar','C','H','H']` | `'ArCH2'` | `'CH2Ar'` (different) |
> | `['Na','C','H','H']` | `'CH2Na'` | `'CH2Na'` (same here) |

---

## CLI

The simplest invocation generates one gas-phase structure with 10 atoms drawn
from elements Z = 1–30:

```bash
pasted --n-atoms 10 --charge 0 --mult 1 \
       --mode gas --region sphere:8 \
       --elements 1-30 --n-samples 1 --seed 42
```

### Placement modes

| Flag | Description | `region` required? |
|---|---|:---:|
| `--mode gas` | Atoms placed uniformly at random inside a sphere or box | ✓ |
| `--mode chain` | Random-walk chain with branching and directional persistence | — |
| `--mode shell` | Central atom surrounded by a coordination shell | — |
| `--mode maxent` | Maximum-entropy placement — neighbors spread as uniformly as possible | ✓ |

> **`region` is required for `gas` and `maxent` modes.**  Use `--region sphere:R`
> (radius R Å) or `--region box:L` (L×L×L Å cube).  The `chain` and `shell`
> modes use their own geometry parameters (`--bond-range`, `--shell-radius`) and
> ignore `region`.

### Generating elongated chain structures

Use `--chain-bias` to bias the chain toward a global axis.
Higher values produce more rod-like structures with larger `shape_aniso`:

```bash
pasted --n-atoms 20 --charge 0 --mult 1 \
       --mode chain --branch-prob 0.0 --chain-bias 0.6 \
       --elements 6,7,8 --n-samples 50 --seed 0 \
       -o chains.xyz
```

### Filtering by disorder metrics

Use `--filter METRIC:MIN:MAX` to keep only structures within a metric range.
Use `-` for an open bound:

```bash
# Keep structures with H_total >= 1.0 and Q6 <= 0.3
pasted --n-atoms 15 --charge 0 --mult 1 \
       --mode gas --region sphere:8 \
       --elements 6,7,8,16 --n-samples 500 --seed 1 \
       --filter "H_total:1.0:-" \
       --filter "Q6:-:0.3" \
       -o filtered.xyz
```

Available metrics: `H_atom`, `H_spatial`, `H_total`, `RDF_dev`,
`shape_aniso`, `Q4`, `Q6`, `Q8`, `graph_lcc`, `graph_cc`,
`ring_fraction`, `charge_frustration`, `moran_I_chi`,
`bond_angle_entropy`, `coordination_variance`, `radial_variance`,
`local_anisotropy`.

```bash
pasted --n-atoms 8 --charge 0 --mult 1 \
       --mode shell --center-z 26 \
       --elements 1-30 --n-samples 20 --seed 7 \
       -o shell_fe.xyz
```

> **Atom count in shell mode:** the number of atoms written to the output file
> may be larger than `--n-atoms` for two reasons.  First, when `--center-z Z`
> is set the center atom is prepended to the coordinate block, adding one atom
> regardless of `n_atoms`.  Second, `add_hydrogen` is enabled by default and
> appends hydrogen atoms to satisfy open valences.  The `n_atoms` parameter
> controls only the number of *shell* atoms placed around the center.  Use
> `--no-add-hydrogen` to suppress hydrogen addition, or inspect
> `Structure.center_sym` to identify the center atom in the output.

### Maximum-entropy mode (`maxent`)

`maxent` places atoms so that their angular distribution around each center is
as uniform as possible.  Like `gas`, it requires a `region` spec:

```bash
pasted --n-atoms 12 --charge 0 --mult 1 \
       --mode maxent --region sphere:6 \
       --elements 6,7,8 --n-samples 50 --seed 42 \
       -o maxent.xyz
```

> **Tip:** `maxent` is slower than `gas` or `chain` because it runs an
> iterative angular-repulsion optimization per structure.  For large N, the
> C++ extension (`HAS_MAXENT_LOOP=True`) is strongly recommended.

For a complete reference of all CLI options, flags, and optimizer mode, see
[CLI reference](cli.md).

---

## Optimizer

`StructureOptimizer` maximizes a user-defined objective function by running
Markov Chain Monte Carlo over atom positions and compositions.  Three methods
are available: Simulated Annealing, Basin-Hopping, and Parallel Tempering.

### Basic usage

```python
from pasted import StructureOptimizer

opt = StructureOptimizer(
    n_atoms=12,
    charge=0,
    mult=1,
    elements="6,7,8,15,16",          # C, N, O, P, S
    objective={"H_total": 1.0, "Q6": -2.0},  # maximize disorder, minimize order
    method="annealing",
    max_steps=5000,
    n_restarts=4,
    seed=42,
)

result = opt.run()
print(result.best)          # best structure found
print(result.summary())     # "restarts=4  best_f=…  worst_f=…  method='annealing'"

for rank, s in enumerate(result, 1):
    print(f"rank {rank}: H_total={s.metrics['H_total']:.3f}  Q6={s.metrics['Q6']:.3f}")
```

### Parallel Tempering

Parallel Tempering runs multiple replicas at different temperatures and
periodically exchanges states, allowing hot replicas to cross energy barriers
while cold replicas refine the best solutions found.

```python
opt = StructureOptimizer(
    n_atoms=12,
    charge=0,
    mult=1,
    elements="6,7,8,15,16",
    objective={"H_total": 1.0, "Q6": -2.0},
    method="parallel_tempering",
    n_replicas=4,               # temperatures: T_end, …, T_start (geometric)
    pt_swap_interval=10,        # attempt replica exchange every 10 steps
    max_steps=2000,
    n_restarts=2,
    T_start=1.0,
    T_end=0.01,
    seed=42,
)

result = opt.run()
print(result.summary())
```

### Electronegativity-targeted optimization

To find structures with frustrated or random electronegativity arrangements:

```python
opt = StructureOptimizer(
    n_atoms=10,
    charge=0,
    mult=1,
    elements="6,7,8,9,14,15,16",     # C, N, O, F, Si, P, S — wide EN range
    objective={
        "charge_frustration": 2.0,    # maximize EN variance across neighbors
        "moran_I_chi": -1.0,          # minimize EN spatial autocorrelation
    },
    method="parallel_tempering",
    n_replicas=4,
    max_steps=3000,
    n_restarts=3,
    seed=7,
)

result = opt.run()
s = result.best
print(f"charge_frustration: {s.metrics['charge_frustration']:.4f}")
print(f"moran_I_chi:        {s.metrics['moran_I_chi']:.3f}")
```

### Affine displacement moves

By default the optimizer uses fragment moves (displacing individual atoms).
Enable **affine moves** to also stretch, compress, and shear the entire
structure — useful for exploring anisotropic configurations:

```python
opt = StructureOptimizer(
    n_atoms=12,
    charge=0,
    mult=1,
    elements="6,7,8",
    objective={"H_total": 1.0},
    allow_affine_moves=True,   # half of displacement moves become affine
    affine_strength=0.15,      # stretch / compress up to ±15 %
    method="annealing",
    max_steps=3000,
    seed=42,
)
result = opt.run()
```

`affine_strength` controls the scale of the transform (default: `0.1`).
Practical range: 0.02–0.4.  Has no effect when `allow_affine_moves=False`.

### Objective function

The objective is **maximized**.  Use negative weights to penalize a metric:

```python
# Dict form: f = sum(w * metric)
objective = {"H_atom": 1.0, "H_spatial": 1.0, "Q6": -2.0}

# Callable form
objective = lambda m: m["H_spatial"] - 2.0 * m["Q6"]
```

Available metrics: all keys in `pasted.ALL_METRICS` —
`H_atom`, `H_spatial`, `H_total`, `RDF_dev`, `shape_aniso`,
`Q4`, `Q6`, `Q8`, `graph_lcc`, `graph_cc`,
`ring_fraction`, `charge_frustration`, `moran_I_chi`,
`bond_angle_entropy`, `coordination_variance`, `radial_variance`,
`local_anisotropy`.

---

### Extended objective with EvalContext: full structure and optimizer state

For objectives that need atomic coordinates, charge/multiplicity, or
optimizer runtime information, use the **2-argument calling convention**.
Any callable with two required positional parameters receives an
[`EvalContext`](api/optimizer.rst) as its second argument:

```python
from pasted import EvalContext, StructureOptimizer
```

The `EvalContext` exposes:

| Field | Type | Description |
|---|---|---|
| `ctx.atoms` | `tuple[str, ...]` | Element symbols |
| `ctx.positions` | `tuple[tuple[float,float,float], ...]` | Coordinates (Å) |
| `ctx.charge` / `ctx.mult` | `int` | Charge and multiplicity |
| `ctx.n_atoms` | `int` | Atom count |
| `ctx.metrics` | `dict[str, float]` | Same as `m` |
| `ctx.to_xyz()` | `str` | XYZ-format string |
| `ctx.step` / `ctx.max_steps` | `int` | MC step index and total |
| `ctx.progress` | `float` | `step / max_steps` ∈ [0, 1) |
| `ctx.temperature` | `float` | Current temperature |
| `ctx.f_current` / `ctx.best_f` | `float` | Current and best scores |
| `ctx.per_atom_q6` | `np.ndarray` | Per-atom Q6, shape `[n_atoms]` |
| `ctx.restart_idx` | `int` | Restart index (0-based) |
| `ctx.element_pool` | `tuple[str, ...]` | Available elements |
| `ctx.cutoff` | `float` | Distance cutoff (Å) |
| `ctx.replica_idx` / `ctx.n_replicas` | `int \| None` | PT-only fields |

**Example 1 — NumPy geometry objective (no external dependencies)**

```python
import numpy as np
from pasted import StructureOptimizer

# 1-arg form: unchanged from before
opt1 = StructureOptimizer(
    n_atoms=10, charge=0, mult=1, elements="6,7,8",
    objective=lambda m: m["H_total"] - 2.0 * m["Q6"],
)

# 2-arg form: use ctx.positions for a purely geometric objective
def rms_spread_objective(m, ctx):
    """Maximize mean pairwise distance."""
    pos = np.array(ctx.positions)
    diffs = pos[:, None, :] - pos[None, :, :]
    dists = np.sqrt((diffs ** 2).sum(axis=-1))
    return float(dists[np.triu_indices(ctx.n_atoms, k=1)].mean())

opt2 = StructureOptimizer(
    n_atoms=10, charge=0, mult=1, elements="6,7,8",
    objective=rms_spread_objective,
    method="annealing", max_steps=2000, seed=42,
)
result = opt2.run()
```

**Example 2 — Adaptive curriculum objective**

```python
def curriculum_objective(m, ctx):
    """Broad exploration early, strong Q6 penalty late."""
    base = m["H_total"]
    if ctx.progress < 0.5:
        return base
    else:
        return base - 3.0 * m["Q6"]

opt = StructureOptimizer(
    n_atoms=15, charge=0, mult=1, elements="6,7,8,16",
    objective=curriculum_objective,
    method="annealing", max_steps=4000, seed=7,
)
```

**Example 3 — Per-atom Q6 locality penalty**

```python
import numpy as np

def local_disorder_objective(m, ctx):
    """Reward Q6 variance; penalize the most locally ordered atom."""
    q6_var = float(np.var(ctx.per_atom_q6))
    q6_max = float(np.max(ctx.per_atom_q6))
    return m["H_total"] + q6_var * 0.5 - q6_max * 1.0

opt = StructureOptimizer(
    n_atoms=20, charge=0, mult=1, elements="6,7,8,15,16",
    objective=local_disorder_objective,
    method="basin_hopping", max_steps=3000, seed=21,
)
```

**Example 4 — xTB single-point energy as objective**

> **External dependency:** `xtb` — GFN-xTB semiempirical tight-binding code (standalone binary).
> Install: `conda install -c conda-forge xtb` (recommended) or `pip install xtb-python`
> Reference: [xtb documentation](https://xtb-docs.readthedocs.io/) | [GitHub](https://github.com/grimme-lab/xtb) | DOI: `10.1021/acs.jctc.8b01176`

```python
import os, subprocess, tempfile

def xtb_energy_objective(m, ctx):
    """Minimize GFN2-xTB single-point energy (negated for maximization).

    Uses ctx.to_xyz() to serialize the structure and ctx.charge / ctx.mult
    for the --chrg / --uhf flags.  Requires the xtb binary on PATH.
    """
    with tempfile.NamedTemporaryFile(
        suffix=".xyz", mode="w", delete=False
    ) as fh:
        fh.write(ctx.to_xyz())
        fname = fh.name
    try:
        proc = subprocess.run(
            ["xtb", fname, "--sp", "--gfn", "2",
             "--chrg", str(ctx.charge),
             "--uhf",  str(ctx.mult - 1)],
            capture_output=True, text=True, timeout=60,
        )
        for line in proc.stdout.splitlines():
            if "TOTAL ENERGY" in line:
                return -float(line.split()[3])   # negate → maximize
        return float("-inf")
    finally:
        os.unlink(fname)

opt = StructureOptimizer(
    n_atoms=8, charge=0, mult=1, elements="6,7,8",
    objective=xtb_energy_objective,
    method="basin_hopping", max_steps=100, seed=7,
)
```

**Example 5 — ASE EMT potential as objective**

> **External dependency:** `ase` — Atomic Simulation Environment.
> Install: `pip install ase`
> Reference: [ASE documentation](https://wiki.fysik.dtu.dk/ase/) | [PyPI](https://pypi.org/project/ase/) | DOI: `10.1088/1361-648X/aa680e`
>
> Note: The EMT calculator below supports only Al, Cu, Ag, Au, Ni, Pd, Pt, H, C, N, O.
> Replace with any ASE-compatible calculator (GPAW, NequIP, MACE, etc.) as needed.

```python
from ase import Atoms
from ase.calculators.emt import EMT

def ase_emt_objective(m, ctx):
    """Use ASE/EMT energy as objective (negated for maximization)."""
    structure = Atoms(
        symbols=list(ctx.atoms),
        positions=list(ctx.positions),
    )
    structure.calc = EMT()
    return -structure.get_potential_energy()   # negate → maximize

opt = StructureOptimizer(
    n_atoms=12, charge=0, mult=1,
    elements="29,79",   # Cu, Au — both supported by EMT
    objective=ase_emt_objective,
    method="annealing", max_steps=500, seed=42,
)
```

---

### `NeighborList` — custom neighbor queries in objectives

`pasted.neighbor_list.NeighborList` is a lightweight, lazily-cached neighbor
list built on `scipy.spatial.cKDTree`.  It is the same data structure that
`compute_all_metrics` uses internally, and it is available as a public API
so that objective functions can compute custom per-bond quantities — such as
Q3 or any other bond-order parameter not in `ALL_METRICS` — without
rebuilding the tree or re-enumerating pairs redundantly.

```python
from pasted.neighbor_list import NeighborList

nl = NeighborList(pts, cutoff)   # build once per evaluation step
```

**Attributes available immediately after construction**

| Attribute | Shape | Description |
|---|---|---|
| `nl.n_atoms` | scalar | Number of atoms N |
| `nl.n_pairs` | scalar | Number of undirected pairs P within cutoff |
| `nl.pairs` | `(P, 2)` | Undirected pair indices, `pairs[k,0] < pairs[k,1]` |
| `nl.rows` / `nl.cols` | `(2P,)` | Directed (bidirectional) source and target indices |
| `nl.dists` | `(P,)` | Undirected pair distances (Å) |
| `nl.all_dists` | `(2P,)` | Directed distances (each undirected distance duplicated) |
| `nl.diff` | `(2P, 3)` | Directed difference vectors `pts[rows] - pts[cols]` |

**Cached properties** (computed on first access, O(1) on subsequent calls)

| Property | Shape | Description |
|---|---|---|
| `nl.deg` | `(N,)` | Per-atom coordination number as float64 |
| `nl.unit_diff` | `(2P, 3)` | Directed unit vectors; coincident pairs (d < 1e-10 Å) use a safe divisor of 1.0 |
| `nl.dists_sq` | `(2P,)` | Squared directed distances `all_dists²` |

**Empty-structure handling:** when `n_atoms=0` or `n_pairs=0` all arrays are
empty and each metric's own early-return guard fires before any arithmetic is
attempted, so no `NaN` or `ZeroDivisionError` is produced.

**Using `NeighborList` inside an objective function**

`ctx.positions` is a `tuple[tuple[float, float, float], ...]` — convert it to
an `ndarray` before passing it to `NeighborList`:

```python
import numpy as np
from pasted.neighbor_list import NeighborList

def my_objective(m, ctx):
    pts = np.array(ctx.positions)          # tuple → ndarray required
    nl = NeighborList(pts, ctx.cutoff)     # built once per step
    # ... use nl.deg, nl.dists, nl.unit_diff, etc.
    return some_score
```

> **Important — set `cutoff` explicitly when using `NeighborList` in
> objectives.**  The default cutoff is `1.5 × median(rᵢ + rⱼ)` over
> covalent radii, which is approximately **2.07 Å** for a C/N/O pool.  At
> this distance only the nearest covalently-bonded neighbors are included
> (mean coordination k ≈ 0.9), which limits the resolution of bond-order
> parameters such as Q3.  For meaningful neighbor statistics, pass an
> explicit cutoff to `StructureOptimizer`:
>
> ```python
> opt = StructureOptimizer(..., cutoff=3.5)   # explicit cutoff in Å
> ```
>
> The same value is then available inside the objective as `ctx.cutoff`, so
> the `NeighborList` and any call to `compute_steinhardt_per_atom` use the
> same distance threshold automatically.

---

### Controlling initial-structure generation retries (`max_init_attempts`)

By default `StructureOptimizer` retries generating the initial structure for
each restart **without a limit** (`max_init_attempts=0`).  This is safe
because the constructor validates at build time that the element pool can
satisfy the charge/multiplicity parity constraint — if construction succeeds,
a valid structure is guaranteed to eventually be found.

Pass a positive integer to cap the number of attempts per restart.  If the
cap is reached without success the restart is skipped and a
`UserWarning` is emitted.  This is useful in automated pipelines where a
wall-time budget matters more than exhaustive search:

```python
from pasted import StructureOptimizer

# At most 100 tries per restart — useful in time-constrained pipelines
opt = StructureOptimizer(
    n_atoms=20,
    charge=0,
    mult=1,
    elements="6,7,8,15,16",
    objective={"H_total": 1.0},
    method="annealing",
    max_steps=3000,
    n_restarts=4,
    max_init_attempts=100,   # cap retries; 0 (default) = unlimited
    seed=42,
)
result = opt.run()
```

Passing an element pool that can *never* satisfy the parity constraint raises
`ValueError` immediately at construction time — no retries are attempted:

```python
# Raises ValueError: all-nitrogen pool cannot satisfy charge=0, mult=1
# N has Z=7 (odd). With an all-odd-Z pool, sum(Z) parity == n_atoms % 2.
# mult=1 requires an even number of electrons, i.e. sum(Z) must be even,
# so n_atoms must be even — but here n_atoms=7 (odd) makes it impossible.
StructureOptimizer(
    n_atoms=7, charge=0, mult=1,
    elements="7",          # nitrogen only — all-odd-Z pool
    objective={"H_total": 1.0},
)
```

---

## GeneratorConfig — immutable configuration object

`GeneratorConfig` is a `frozen=True` dataclass that encapsulates every
parameter of `StructureGenerator`.  It gives full mypy / IDE type-checking
and allows safe one-field overrides via `dataclasses.replace`.

```python
import dataclasses
from pasted import GeneratorConfig, StructureGenerator

# Build once
cfg = GeneratorConfig(
    n_atoms=20, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8", n_samples=100, seed=42,
)

# Pass to the class API
result = StructureGenerator(cfg).generate()

# One-field override — creates a new config, does not mutate the original
cfg_new_seed = dataclasses.replace(cfg, seed=99)
result2 = StructureGenerator(cfg_new_seed).generate()
```

`generate()` also accepts a config directly:

```python
from pasted import generate, GeneratorConfig

result = generate(GeneratorConfig(n_atoms=12, charge=0, mult=1,
                                  mode="chain", elements="6,7,8",
                                  n_samples=50, seed=0))
```

The original keyword-argument style (`StructureGenerator(n_atoms=..., ...)`
and `generate(n_atoms=..., ...)`) is **fully backward-compatible** and
continues to work unchanged.  All instance attributes
(`gen.n_atoms`, `gen.seed`, …) are still accessible via `__getattr__` proxy.

---

## Affine transforms in StructureGenerator

Setting `affine_strength > 0` applies a random affine transformation —
stretch/compress along one axis, shear one axis pair, and optionally a
small per-atom jitter — to each generated structure **before**
`relax_positions`.  The center of mass is pinned throughout, so the
structure stays centered.  The transform creates more anisotropic initial
geometries while still guaranteeing clash-free output after relaxation.

Three sub-operations can be controlled independently via
`affine_stretch`, `affine_shear`, and `affine_jitter`.  When any of
these is `None` (the default), it falls back to `affine_strength`.
Setting one of them to `0.0` disables that operation entirely regardless
of `affine_strength`.

### Basic usage

```python
from pasted import generate

# Global control: all three operations use the same strength
result = generate(
    n_atoms=20, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8", n_samples=50, seed=42,
    affine_strength=0.2,   # ±20 % stretch, ±10 % shear, jitter ∝ 0.2
)
```

### Per-operation control

```python
from pasted import generate

# Strong axial elongation only — no shear, no per-atom noise
result_stretch = generate(
    n_atoms=20, charge=0, mult=1,
    mode="chain", elements="6,7,8",
    n_samples=50, seed=0,
    affine_strength=0.2,
    affine_stretch=0.4,   # large stretch/compress (±40 %)
    affine_shear=0.0,     # shear disabled
    affine_jitter=0.0,    # per-atom jitter disabled
)

# Pure shear distortion — no axial scaling
result_shear = generate(
    n_atoms=20, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8", n_samples=50, seed=0,
    affine_strength=0.2,
    affine_stretch=0.0,   # stretch disabled
    affine_shear=0.3,     # shear only
    affine_jitter=0.0,
)

# Per-atom jitter only — break symmetry without global distortion
result_jitter = generate(
    n_atoms=20, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8", n_samples=50, seed=0,
    affine_strength=0.1,
    affine_stretch=0.0,
    affine_shear=0.0,
    affine_jitter=0.3,    # fine-grain per-atom noise only
)
```

### Parameter reference

| Parameter | Controls | Range | Default |
|---|---|---|---|
| `affine_strength` | Global fallback for all three operations | 0.0–0.4 | `0.0` (disabled) |
| `affine_stretch` | Scale factor along one random axis: `Uniform(1−s, 1+s)` | 0.0–0.4 | `None` (→ `affine_strength`) |
| `affine_shear` | Off-diagonal matrix element: `Uniform(-s/2, s/2)` | 0.0–0.4 | `None` (→ `affine_strength`) |
| `affine_jitter` | Per-atom translation noise ∝ `move_step × s` | 0.0–0.4 | `None` (→ `affine_strength`) |

> **Note on `affine_jitter` in StructureGenerator:** when used inside
> `generate()` or `StructureGenerator`, the internal `move_step` is `0.0`,
> so `affine_jitter` has no visible effect — the jitter term is gated on
> `move_step > 0`.  `affine_jitter` is therefore only meaningful in
> `StructureOptimizer` (where `move_step` is set per MC step).

Practical `affine_strength` guide:

| Value | Effect |
|---|---|
| `0.0` (default) | disabled — backward-compatible, no transform |
| `0.05`–`0.1` | subtle anisotropy, negligible overhead |
| `0.2`–`0.3` | noticeable elongation or shear |
| `0.4` | strong distortion; may require larger `region` to avoid clashes |

Works across **all placement modes** (`gas`, `chain`, `shell`, `maxent`).
CLI flags: `--affine-strength S`, `--affine-stretch S`, `--affine-shear S`,
`--affine-jitter S`.

```bash
# CLI: axial elongation only, no shear
pasted --n-atoms 20 --charge 0 --mult 1 \
       --mode chain --elements 6,7,8 --n-samples 50 --seed 0 \
       --affine-strength 0.2 --affine-stretch 0.4 \
       --affine-shear 0.0 --affine-jitter 0.0
```

### Affine moves in StructureOptimizer

`StructureOptimizer` supports the same four parameters through
`allow_affine_moves=True`.  When enabled, half of all displacement moves
are replaced by affine moves, allowing the optimizer to explore
elongated, compressed, or sheared configurations that fragment moves
cannot reach efficiently.

```python
from pasted import StructureOptimizer

# Strong axial stretch — useful for maximizing shape_aniso
opt = StructureOptimizer(
    n_atoms=20, charge=0, mult=1,
    elements="6,7,8",
    objective={"shape_aniso": 2.0, "H_total": 1.0},
    allow_affine_moves=True,
    affine_strength=0.2,
    affine_stretch=0.4,   # dominant stretch moves
    affine_shear=0.0,     # no shear
    affine_jitter=0.0,    # no per-atom noise
    method="annealing",
    max_steps=2000, n_restarts=2, seed=0,
)
result = opt.run()
print(f"shape_aniso = {result.best.metrics['shape_aniso']:.4f}")

# Combined stretch + shear for complex anisotropy
opt2 = StructureOptimizer(
    n_atoms=20, charge=0, mult=1,
    elements="6,7,8",
    objective={"shape_aniso": 2.0, "H_total": 1.0},
    allow_affine_moves=True,
    affine_strength=0.2,
    affine_stretch=0.25,
    affine_shear=0.15,
    affine_jitter=0.0,    # disable noise for pure geometric distortion
    method="annealing",
    max_steps=2000, n_restarts=2, seed=0,
)
result2 = opt2.run()
print(f"shape_aniso = {result2.best.metrics['shape_aniso']:.4f}")
```

> **Relationship to `StructureGenerator`:** both classes use the same
> `_affine_move` function from `pasted._placement`.  In `StructureGenerator`
> the transform is applied once per structure before `relax_positions`; in
> `StructureOptimizer` it is applied at each accepted MC step when
> `allow_affine_moves=True`.  The `affine_jitter` parameter only has a
> visible effect in the Optimizer (where `move_step > 0`).

---

## Element sampling control (StructureGenerator)

### Element pool specification

The `elements=` parameter (Python API) and `--elements` (CLI) control which
elements are available for random sampling.  Three forms are supported:

| Form | Example | Meaning |
|---|---|---|
| Atomic-number spec string | `"6,7,8"` | C, N, O (by atomic number Z) |
| Atomic-number range string | `"1-30"` | H through Zn |
| Combined ranges + singles | `"1-10,26,28"` | H–Ne plus Fe and Ni |
| List of element symbols | `["C", "N", "O"]` | explicit symbol list |
| `None` (omitted) | — | all Z = 1–106 |

> **Important:** when passing a *string*, it must contain **atomic numbers
> (integers)**, not element symbols.  `elements="C,N,O"` raises `ValueError`;
> use `elements="6,7,8"` or `elements=["C", "N", "O"]` instead.

> **Note (v0.3.5):** `parse_element_spec()` now accepts a `list[str]` of
> element symbols directly (e.g. `["C", "N", "O"]`).  Previously, calling
> `parse_element_spec(["C", "N", "O"])` raised `AttributeError` because
> the function attempted to call `.split(",")` on the list.  The fix has
> no effect on the `StructureGenerator` / `generate()` keyword API, which
> already handled symbol lists correctly via an internal branch.

```python
# Correct — numeric atomic-number string
gen = StructureGenerator(n_atoms=10, charge=0, mult=1, mode="chain",
                         elements="6,7,8")

# Correct — explicit symbol list (passed to StructureGenerator or parse_element_spec)
gen = StructureGenerator(n_atoms=10, charge=0, mult=1, mode="chain",
                         elements=["C", "N", "O"])

from pasted._atoms import parse_element_spec
assert parse_element_spec(["C", "N", "O"]) == ["C", "N", "O"]  # now works

# WRONG — symbol string raises ValueError
# gen = StructureGenerator(..., elements="C,N,O")  # ← ValueError!
```

### Biased element fractions

By default each element is sampled with equal probability.  Pass
`element_fractions` to shift the distribution:

```python
from pasted import StructureGenerator

gen = StructureGenerator(
    n_atoms=20, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8",
    element_fractions={"C": 0.6, "N": 0.3, "O": 0.1},  # C-rich
    n_samples=50, seed=0,
)
result = gen.generate()
```

Weights are *relative* — they are normalized internally.
`{"C": 6, "N": 3, "O": 1}` is equivalent to the above.
Elements absent from the dict (but present in the pool) receive weight `1.0`.

> **Note — keys must be within the element pool:** every key in
> `element_fractions` must correspond to an element present in the `elements=`
> pool.  Passing a symbol that is not in the pool raises `ValueError` at
> construction time.  This is by design: a weight for an element that can
> never be sampled would silently mislead callers about the actual sampling
> distribution.

### Element count bounds

Use `element_min_counts` and `element_max_counts` to guarantee or cap the
number of specific atoms in each generated structure:

```python
gen = StructureGenerator(
    n_atoms=15, charge=0, mult=1,
    mode="gas", region="sphere:10",
    elements="6,7,8,15,16",          # C, N, O, P, S
    element_min_counts={"C": 4},      # always at least 4 carbons
    element_max_counts={"N": 3, "O": 3},  # cap nitrogen and oxygen
    n_samples=100, seed=42,
)
result = gen.generate()

from collections import Counter
for s in result:
    c = Counter(s.atoms)
    assert c["C"] >= 4
    assert c.get("N", 0) <= 3
    assert c.get("O", 0) <= 3
```

Bounds are enforced during atom sampling, before the parity check.
A `ValueError` is raised at construction time when constraints are
inconsistent (sum of mins > `n_atoms`, or min > max for any element).

### Combining fractions and bounds

All three parameters can be used together:

```python
gen = StructureGenerator(
    n_atoms=12, charge=0, mult=1,
    mode="chain",
    elements="6,7,8",
    element_fractions={"C": 5, "N": 2, "O": 1},
    element_min_counts={"C": 2},
    element_max_counts={"N": 4},
    n_samples=30, seed=7,
)
```

---

## Position-only optimization (StructureOptimizer)

Set `allow_composition_moves=False` to fix the composition and only
optimize atomic positions.  This is useful when the stoichiometry is
predetermined:

```python
from pasted import StructureOptimizer, Structure

# Load a structure with a fixed composition
initial = Structure.from_xyz("my_structure.xyz")

opt = StructureOptimizer(
    n_atoms=len(initial),
    charge=initial.charge,
    mult=initial.mult,
    elements=list(set(initial.atoms)),
    objective={"H_total": 1.0, "Q6": -2.0},
    allow_composition_moves=False,   # position-only
    method="annealing",
    max_steps=5000,
    seed=42,
)

result = opt.run(initial=initial)
print(result.best)
# Composition is identical to initial; only positions have changed
assert sorted(result.best.atoms) == sorted(initial.atoms)
```

---

## Composition-only optimization (StructureOptimizer)

Set `allow_displacements=False` to fix the atomic coordinates and only
optimize element types.  This is useful when exploring compositional
disorder on a pre-relaxed geometry (e.g. a fixed lattice).

The element pool does not need to match the initial structure's composition.
Any atoms outside the pool are automatically replaced by parity-compatible
pool elements before the MC loop begins, so cross-pool optimization works
reliably with all three methods (`"annealing"`, `"basin_hopping"`, and
`"parallel_tempering"`):

```python
from pasted import StructureOptimizer, generate

# Starting geometry can use any element set — the optimizer will sanitize
# foreign atoms to the target pool before the first MC step.
initial = generate(
    n_atoms=10, charge=0, mult=1,
    mode="gas", region="sphere:8",
    elements="6,7,8",          # C / N / O geometry
    n_samples=50, seed=0,
)[0]

opt = StructureOptimizer(
    n_atoms=len(initial),
    charge=initial.charge,
    mult=initial.mult,
    elements=["Cr", "Mn", "Fe", "Co", "Ni"],  # Cantor alloy pool
    objective={
        "moran_I_chi": -1.0,          # minimize EN spatial autocorrelation
        "charge_frustration": 2.0,    # maximize EN variance across neighbors
    },
    allow_displacements=False,   # composition-only; coordinates fixed
    method="annealing",
    max_steps=5000,
    seed=42,
)

result = opt.run(initial=initial)
print(result.best.comp)   # e.g. 'CoCr3Fe3MnNi2'
# Positions are identical to initial; only element labels have changed
import numpy as np
np.testing.assert_allclose(
    np.array(result.best.positions), np.array(initial.positions)
)
```

> **Choosing an objective for composition-only runs:** each composition move
> selects a random atom and replaces it with a different element drawn from
> the pool, which changes the element-count histogram.  However, metrics
> that capture the *spatial arrangement* of elements — `moran_I_chi`
> (electronegativity autocorrelation) and `charge_frustration` (EN variance
> across neighbor pairs) — respond most sensitively to which element sits
> where and are therefore the recommended choices for composition-only
> optimization.  Pure composition-count metrics such as `H_atom` tend to
> plateau quickly once the pool coverage is broad.

> **Note**: `allow_displacements=False` and `allow_composition_moves=False`
> cannot both be set — at least one move type must be enabled.  Attempting
> to do so raises a `ValueError`.

---

## Optimizer case studies

### Case study 1 — Reproducing a target disorder profile

This example shows how to drive `StructureOptimizer` toward a reference
structure defined by a known set of disorder metrics.  The target is a
dense 40-atom C/N sphere with high Moran's I spatial autocorrelation
(`moran_I_chi ≈ 0.94`), a fully connected graph (`graph_lcc = 1.0`), and
a high ring fraction (`ring_fraction = 0.875`).  The `cutoff` parameter is
set to match the value used when the reference metrics were computed.

```python
from pasted import StructureOptimizer

# Target metrics (from moran1_dense_sphere.xyz, computed at cutoff=5.0 Å):
#   moran_I_chi   = 0.9446   graph_lcc     = 1.0000
#   ring_fraction = 0.8750   charge_frustration = 0.0058
#   H_spatial     = 3.4420   Q6            = 0.3335

opt = StructureOptimizer(
    n_atoms=40,
    charge=0,
    mult=1,
    elements=["C", "N"],
    objective={
        "moran_I_chi":        3.0,   # primary target: high spatial EN autocorrelation
        "ring_fraction":      2.0,   # dense ring connectivity
        "graph_lcc":          2.0,   # fully connected graph
        "H_spatial":          1.0,   # spatial disorder
        "charge_frustration": -5.0,  # suppress (target ≈ 0)
    },
    cutoff=5.0,                      # match the reference metric cutoff
    method="parallel_tempering",
    n_replicas=4,
    pt_swap_interval=10,
    max_steps=3000,
    n_restarts=3,
    T_start=1.0,
    T_end=0.01,
    seed=42,
)

result = opt.run()
best = result.best
print(result.summary())
print(f"moran_I_chi   = {best.metrics['moran_I_chi']:.4f}  (target 0.9446)")
print(f"ring_fraction = {best.metrics['ring_fraction']:.4f}  (target 0.8750)")
print(f"graph_lcc     = {best.metrics['graph_lcc']:.4f}  (target 1.0000)")
best.write_xyz("similar_to_target.xyz")
```

**Tips for target-matching runs:**

- Always pass `cutoff=` explicitly so that PASTED uses the same distance
  threshold as was used to compute the reference metrics.
- Use `n_restarts ≥ 3` with `method="parallel_tempering"` — the landscape
  is rugged and hot replicas help escape local optima.
- Negative weights (`charge_frustration: -5.0`) suppress metrics that
  should stay near zero; large magnitude is important to counteract the
  positive objectives.

---

### Case study 2 — Geometry search with an ASE calculator

`EvalContext` lets you call any ASE-compatible calculator inside the
objective function, connecting PASTED's stochastic search to the entire
ASE ecosystem (xTB, MACE, NequIP, GPAW, …).  This example uses the
built-in EMT potential to search for a low-energy methane-like (CH₄)
geometry.

> **External dependency:** `ase` — Atomic Simulation Environment.
> Install: `pip install ase`
>
> **Note:** the EMT calculator supports only Al, Cu, Ag, Au, Ni, Pd, Pt,
> H, C, N, O.  Replace with any ASE-compatible calculator as needed.

```python
from ase import Atoms
from ase.calculators.emt import EMT
from ase.optimize import BFGS

from pasted import EvalContext, StructureGenerator, StructureOptimizer

# ── Step 1: generate C₁H₄ initial geometries with PASTED ────────────────
initial_structs = StructureGenerator(
    n_atoms=5,
    charge=0,
    mult=1,
    mode="gas",
    region="sphere:3",
    elements=["C", "H"],
    element_min_counts={"C": 1},
    element_max_counts={"C": 1},
    n_samples=100,
    seed=0,
).generate()

initial = initial_structs[0]   # pick the first C₁H₄ structure

# ── Step 2: basin-hopping with ASE/EMT objective ─────────────────────────
def ase_emt_objective(m: dict, ctx: EvalContext) -> float:
    """Maximize stability — lower EMT potential energy is better."""
    structure = Atoms(
        symbols=list(ctx.atoms),
        positions=list(ctx.positions),
    )
    structure.calc = EMT()
    try:
        return -structure.get_potential_energy()   # negate: maximize = minimize energy
    except Exception:
        return float("-inf")

opt = StructureOptimizer(
    n_atoms=len(initial),
    charge=initial.charge,
    mult=initial.mult,
    elements=list(set(initial.atoms)),
    objective=ase_emt_objective,
    allow_composition_moves=False,   # fix C₁H₄ stoichiometry
    method="basin_hopping",
    max_steps=500,
    n_restarts=3,
    seed=7,
)
result = opt.run(initial=initial)
best = result.best
print(result.summary())

# ── Step 3: final geometry relaxation with ASE BFGS ──────────────────────
ase_mol = Atoms(
    symbols=list(best.atoms),
    positions=list(best.positions),
)
ase_mol.calc = EMT()
bfgs = BFGS(ase_mol, logfile=None)
bfgs.run(fmax=0.01)

print(f"EMT energy after BFGS: {ase_mol.get_potential_energy():.4f} eV")
for sym, pos in zip(ase_mol.get_chemical_symbols(), ase_mol.positions):
    print(f"  {sym:2s}  {pos[0]:8.4f}  {pos[1]:8.4f}  {pos[2]:8.4f}")
```

**Workflow summary:**

1. `StructureGenerator` samples random clash-free starting geometries
   with the required stoichiometry.
2. `StructureOptimizer` explores configuration space using PASTED's MC
   moves, calling the ASE calculator at each accepted step.
3. The best PASTED geometry is passed to ASE `BFGS` for a final
   gradient-based refinement.

This three-stage pipeline (sample → stochastic search → local minimize)
is a general pattern for any external potential.

---

### Case study 3 — Two-phase curriculum objective

`EvalContext.progress` returns the fractional progress of the current
restart (range `[0.0, 1.0)`), enabling objectives that change behavior
over the course of a run.  This is useful for annealing-style curricula
where early exploration and late refinement require different signals.

```python
from pasted import EvalContext, StructureOptimizer

def curriculum_objective(m: dict, ctx: EvalContext) -> float:
    """
    Two-phase curriculum driven by optimizer progress:

    Phase 1 (first 50% of steps):
        Maximize spatial disorder only — broad exploration of
        configuration space without penalizing ordered regions.

    Phase 2 (last 50% of steps):
        Continue maximizing spatial disorder AND actively suppress
        crystalline order (Q6).  Switching late avoids trapping
        the search in disordered-but-ordered local minima early on.
    """
    if ctx.progress < 0.5:
        return m["H_spatial"] * 2.0
    else:
        return m["H_spatial"] * 2.0 - m["Q6"] * 3.0

opt = StructureOptimizer(
    n_atoms=20,
    charge=0,
    mult=1,
    elements="6,7,8,15,16",          # C, N, O, P, S
    objective=curriculum_objective,
    method="annealing",
    max_steps=2000,
    n_restarts=4,
    seed=99,
)

result = opt.run()
best = result.best
print(result.summary())
print(f"H_spatial = {best.metrics['H_spatial']:.4f}")
print(f"Q6        = {best.metrics['Q6']:.4f}")
```

**When to use curriculum objectives:**

- When the search landscape has a conflict between early exploration and
  late refinement (e.g., maximizing disorder while also minimizing a
  specific metric).
- When you want to prevent the optimizer from committing too early to a
  particular region of metric space.
- Use `ctx.step` for step-based schedules and `ctx.temperature` to tie
  the curriculum to the annealing schedule directly.

---

### Case study 4 — Custom bond-order metric (Q3) via `NeighborList`

`ALL_METRICS` provides Q4, Q6, and Q8, but Steinhardt parameters for other
angular momenta — such as Q3, which is sensitive to three-fold local
symmetry — can be computed inside a 2-argument objective using
`NeighborList` and `compute_steinhardt_per_atom`.  This pattern generalises
to any bond-level quantity: just build a `NeighborList` from `ctx.positions`
and operate on its arrays directly.

**Why Q3?**  Q3 is nonzero for structures with three-fold (triangular or
kagomé-like) local coordination.  It is absent from `ALL_METRICS` because
it is rarely needed for the fuzzing workflows PASTED targets, but it is
straightforward to add via the objective callable.

```python
import numpy as np
from pasted import EvalContext, StructureOptimizer
from pasted.neighbor_list import NeighborList
from pasted._metrics import compute_steinhardt_per_atom

# ── Objective: maximize mean Q3 across all atoms ─────────────────────────
def q3_objective(m: dict, ctx: EvalContext) -> float:
    """Maximize the Steinhardt Q3 bond-order parameter.

    Uses NeighborList to avoid rebuilding the cKDTree independently from
    the per-atom Steinhardt call.  ctx.cutoff is supplied by the optimizer
    and is the same cutoff used for all built-in metrics, ensuring consistency.
    """
    pts = np.array(ctx.positions)                                   # tuple → ndarray
    q3_per_atom = compute_steinhardt_per_atom(pts, [3], ctx.cutoff)["Q3"]
    return float(q3_per_atom.mean())

# ── Combined objective: maximize Q3 and compositional disorder ───────────
def q3_disorder_objective(m: dict, ctx: EvalContext) -> float:
    """Maximize Q3 and spatial disorder simultaneously."""
    pts = np.array(ctx.positions)
    q3_per_atom = compute_steinhardt_per_atom(pts, [3], ctx.cutoff)["Q3"]
    return float(q3_per_atom.mean()) + m["H_total"]

# ── Run the optimizer ─────────────────────────────────────────────────────
# Pass cutoff=3.5 explicitly so that ctx.cutoff gives k_avg ≈ 1.5 rather
# than the default ~2.07 Å (k_avg ≈ 0.9), which is too sparse for Q3.
opt = StructureOptimizer(
    n_atoms=20,
    charge=0,
    mult=1,
    elements="6,7,8",
    objective=q3_disorder_objective,
    cutoff=3.5,           # explicit cutoff — essential for meaningful Q3
    method="annealing",
    max_steps=2000,
    n_restarts=3,
    seed=42,
)

result = opt.run()
best = result.best
print(result.summary())

# Verify Q3 independently at the same cutoff
q3_check = compute_steinhardt_per_atom(
    np.array(best.positions), [3], 3.5
)["Q3"]
print(f"Q3 (global mean) = {float(q3_check.mean()):.4f}")
print(f"H_total          = {best.metrics['H_total']:.4f}")
print(f"comp             = {best.comp}")
```

**Extending to other `NeighborList`-based quantities**

The same pattern works for any bond-level descriptor.  For example, use
`nl.unit_diff` and `nl.deg` to build a per-atom angular anisotropy score,
or `nl.dists` to compute a custom radial moment:

```python
def custom_radial_objective(m: dict, ctx: EvalContext) -> float:
    """Penalize structures where neighbor distances are too uniform."""
    pts = np.array(ctx.positions)
    nl = NeighborList(pts, ctx.cutoff)
    if nl.n_pairs == 0:
        return 0.0
    # per-atom variance of neighbor distances
    per_atom_var = np.zeros(nl.n_atoms)
    for i in range(nl.n_atoms):
        mask = nl.rows == i
        if mask.sum() > 1:
            per_atom_var[i] = float(np.var(nl.all_dists[mask]))
    return float(per_atom_var.mean()) + m["H_total"]
```

**Tips for custom-metric objectives**

- Always pass `cutoff=` to `StructureOptimizer` and use `ctx.cutoff` inside
  the objective so that the neighbor list and the built-in metrics share the
  same distance threshold.
- Guard against `nl.n_pairs == 0` (all atoms beyond cutoff) to avoid
  division by zero or empty-array errors; returning `0.0` or `float("-inf")`
  is appropriate depending on whether zero connectivity should be penalized.
- `NeighborList.deg`, `unit_diff`, and `dists_sq` are cached on first access
  within a single `NeighborList` instance, so calling them multiple times in
  one objective evaluation is free.