diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index a0bc793..fc2229d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11"]
+        python-version: ["3.10", "3.11", "3.12"]
     steps:
       - uses: actions/checkout@v4
 
@@ -27,11 +27,29 @@ jobs:
           # Lightweight stack covering the MSM and CV-method tests. The heavy MD
           # backends (openmm, MDAnalysis) are not required: tests that need them
           # skip automatically via pytest.importorskip.
-          pip install numpy scipy scikit-learn pydantic pyyaml deeptime pytest ruff
+          pip install numpy scipy scikit-learn shapely pydantic pyyaml deeptime pytest ruff
           pip install torch --index-url https://download.pytorch.org/whl/cpu
 
-      - name: Lint (ruff)
-        run: ruff check autosampler/msm autosampler/execution autosampler/analysis autosampler/analysis_cli.py autosampler/init_cli.py autosampler/templates.py autosampler/binning/we.py autosampler/binning/adaptive.py autosampler/spaces/registry.py autosampler/spaces/spib.py autosampler/spaces/feature_selection.py autosampler/spaces/retraining.py autosampler/utils/seeds.py autosampler/spawners/we.py autosampler/spawners/msm.py autosampler/reporting.py autosampler/engines/base.py tests
+      - name: Lint (ruff, whole tree)
+        run: ruff check autosampler tests
 
       - name: Run test suite
         run: pytest -q
+
+  docs:
+    name: Docs build
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+      - name: Install docs dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install mkdocs-material "mkdocstrings[python]"
+          # Make the package importable for mkdocstrings/griffe without pulling
+          # the heavy runtime deps (griffe analyses the source statically).
+          pip install -e . --no-deps
+      - name: Build site
+        run: mkdocs build
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
new file mode 100644
index 0000000..a0d0876
--- /dev/null
+++ b/.github/workflows/release.yml
@@ -0,0 +1,46 @@
+name: Release
+
+# Publishes to PyPI when a version tag (v*) is pushed.
+#
+# Uses PyPI Trusted Publishing (OIDC) — no API token/secret required. Configure a
+# trusted publisher for this repository + workflow on PyPI before the first
+# release: https://docs.pypi.org/trusted-publishers/
+on:
+  push:
+    tags: ["v*"]
+
+jobs:
+  build:
+    name: Build distributions
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+      - name: Build sdist and wheel
+        run: |
+          python -m pip install --upgrade pip build
+          python -m build
+      - name: Check metadata
+        run: |
+          pip install twine
+          twine check dist/*
+      - uses: actions/upload-artifact@v4
+        with:
+          name: dist
+          path: dist/
+
+  publish:
+    name: Publish to PyPI
+    needs: build
+    runs-on: ubuntu-latest
+    environment: pypi
+    permissions:
+      id-token: write  # required for trusted publishing
+    steps:
+      - uses: actions/download-artifact@v4
+        with:
+          name: dist
+          path: dist/
+      - uses: pypa/gh-action-pypi-publish@release/v1
diff --git a/AutoSampler_Review.pdf b/AutoSampler_Review.pdf
new file mode 100644
index 0000000..472b8d6
Binary files /dev/null and b/AutoSampler_Review.pdf differ
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fd56b2f..fffe02e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -122,6 +122,26 @@ adds a flux-weighted **transition-matrix convergence** gate with
 - Made `autosampler.spaces` import lazily so lightweight modules (e.g. the CV
   registry) import without MDAnalysis / torch.
 
+### Production-readiness hardening (pre-release pass)
+- **Reproducibility:** threaded the configured seed into all learned-CV training
+  (SPIB no longer hardcodes `seed=42`; the torch RNG is reseeded before every
+  `fit`), and `SeedManager` now requests deterministic torch algorithms.
+- **Robustness:** the local backend tolerates a single walker's failure (and adds
+  an opt-in `execution.walker_timeout` hang guard) instead of aborting the batch;
+  delta-checkpoint resume reconstructs the full history (fixing a truncated
+  `autosampler-path`), writes atomically, and tolerates a corrupt delta; fixed a
+  target-mode spawn crash and a deep-TICA device mismatch.
+- **Packaging:** OpenMM is now an optional, lazily-imported backend so the base
+  `pip install` resolves; full PyPI metadata + single-sourced version;
+  `CITATION.cff`; tag-driven PyPI release workflow.
+- **CI/quality:** lint the whole tree (was a hand-picked subset), test on Python
+  3.10–3.12, build the docs in CI.
+- **Docs/examples/tests:** a self-contained CPU-only alanine-dipeptide
+  hello-world; example configs for SPIB / deep-TICA / WE / target / PBS; an
+  examples index; an API reference, references/citations page, and full CLI
+  reference; +23 tests (delta checkpoint, reproducibility, timeout, spawners,
+  paths). Suite 95 → 118.
+
 ### Fixed (Phase 2)
 - Renamed `AdaptiveSpaceModel.fited` → `fitted` (with a backwards-compatible
   loader for old checkpoints).
diff --git a/CITATION.cff b/CITATION.cff
new file mode 100644
index 0000000..d9908dc
--- /dev/null
+++ b/CITATION.cff
@@ -0,0 +1,36 @@
+cff-version: 1.2.0
+message: "If you use AutoSampler in your research, please cite it as below."
+title: "AutoSampler: an MSM-convergence-driven adaptive sampling framework for molecular dynamics"
+abstract: >-
+  AutoSampler is a modular framework for autonomous adaptive molecular-dynamics
+  sampling. It runs short MD walkers, projects frames into fixed or machine-learned
+  collective-variable spaces, restarts walkers from informative regions, builds a
+  Markov State Model each iteration, and stops when the MSM has converged.
+type: software
+# TODO(release): confirm the full author list, ORCIDs, and affiliations before
+# submitting the publication. The entries below are drawn from the repository and
+# must be verified.
+authors:
+  - family-names: Maity
+    given-names: Dibyendu
+    email: dibyendumaity1999@bose.res.in
+    # orcid: "https://orcid.org/0000-0000-0000-0000"
+    # affiliation: "S. N. Bose National Centre for Basic Sciences"
+  - family-names: Chakrabarty
+    given-names: Suman
+    email: chakrabarty.suman@gmail.com
+    # orcid: "https://orcid.org/0000-0000-0000-0000"
+    # affiliation: "S. N. Bose National Centre for Basic Sciences"
+version: 2.0.0
+license: MIT
+repository-code: "https://github.com/TeamSuman/AutoSampler"
+url: "https://github.com/TeamSuman/AutoSampler"
+keywords:
+  - molecular dynamics
+  - adaptive sampling
+  - Markov state model
+  - collective variables
+  - enhanced sampling
+# TODO(release): add the DOI once a Zenodo archive / preprint is minted, and add a
+# `preferred-citation:` block pointing to the journal article when published.
+# doi: "10.5281/zenodo.XXXXXXX"
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..a2b40d1
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,3 @@
+include README.md CHANGELOG.md LICENSE CITATION.cff env.yml
+recursive-include docs *.md
+exclude AutoSampler_Changelog.pdf AutoSampler_Review.pdf
diff --git a/README.md b/README.md
index fa0b6e1..1ebc626 100644
--- a/README.md
+++ b/README.md
@@ -350,3 +350,9 @@ This is an active research codebase. The examples are useful starting points,
 but production scientific claims should be made only after checking CV
 definitions, sampling bounds, lineage connectivity, and system-specific
 validation.
+
+## How to cite
+
+If you use AutoSampler in your research, please cite it. Citation metadata is in
+[`CITATION.cff`](CITATION.cff) (GitHub renders a "Cite this repository" button
+from it). A DOI and the accompanying publication will be added on release.
diff --git a/autosampler/binning/spatial.py b/autosampler/binning/spatial.py
index 7d2eec2..5f1dfc0 100644
--- a/autosampler/binning/spatial.py
+++ b/autosampler/binning/spatial.py
@@ -372,7 +372,7 @@ def _finite_voronoi_regions_2d(
         radius = float(np.ptp(vor.points, axis=0).max() * 2.0)
 
     all_ridges: dict[int, list[tuple[int, int, int]]] = {}
-    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
+    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices, strict=False):
         all_ridges.setdefault(p1, []).append((p2, v1, v2))
         all_ridges.setdefault(p2, []).append((p1, v1, v2))
 
@@ -405,7 +405,7 @@ def _finite_voronoi_regions_2d(
         vs = np.asarray([new_vertices[vertex] for vertex in new_region])
         centroid = vs.mean(axis=0)
         angles = np.arctan2(vs[:, 1] - centroid[1], vs[:, 0] - centroid[0])
-        new_region = [vertex for _, vertex in sorted(zip(angles, new_region))]
+        new_region = [vertex for _, vertex in sorted(zip(angles, new_region, strict=False))]
         new_regions.append(new_region)
 
     return new_regions, np.asarray(new_vertices, dtype=float)
diff --git a/autosampler/checkpoints/manager.py b/autosampler/checkpoints/manager.py
index 7d74e77..23b996e 100644
--- a/autosampler/checkpoints/manager.py
+++ b/autosampler/checkpoints/manager.py
@@ -1,8 +1,10 @@
 import logging
+import os
 import pickle
-import torch
 from pathlib import Path
-from typing import Any, Dict, Tuple
+from typing import Any
+
+import torch
 
 # On-disk checkpoint format version. Bump when the layout changes; ``load``
 # tolerates older checkpoints (a missing version file is treated as v1).
@@ -13,6 +15,50 @@
 _TORCH_ENCODER_MODES = ("tvae", "vampnet", "spib")
 
 
+def _atomic_pickle(obj: Any, path: Path) -> None:
+    """Pickle ``obj`` to ``path`` atomically (write tmp, then os.replace).
+
+    Prevents a crash mid-write (e.g. an HPC walltime kill) from leaving a
+    truncated file — important for delta history, where every checkpoint's
+    ``history.pkl`` is needed to reconstruct the full history on resume.
+    """
+    tmp = path.with_name(path.name + ".tmp")
+    with open(tmp, "wb") as handle:
+        pickle.dump(obj, handle)
+    os.replace(tmp, path)
+
+
+def reconstruct_history(checkpoint_root: Path, iteration: int) -> dict[Any, Any]:
+    """Merge the per-checkpoint delta ``history.pkl`` files (for all iterations
+    ``<= iteration``) into the full cumulative history.
+
+    Each key normally lives in exactly one delta; on overlap the newer checkpoint
+    wins. Unreadable deltas (truncated by a crash) are skipped with a warning
+    rather than aborting the whole restore.
+    """
+    iters = sorted(
+        int(path.name.removeprefix("iter_"))
+        for path in checkpoint_root.glob("iter_*")
+        if path.is_dir()
+        and path.name.removeprefix("iter_").isdigit()
+        and int(path.name.removeprefix("iter_")) <= iteration
+    )
+    full: dict[Any, Any] = {}
+    for it in iters:
+        hist_file = checkpoint_root / f"iter_{it}" / "history.pkl"
+        if not hist_file.exists():
+            continue
+        try:
+            with open(hist_file, "rb") as handle:
+                part = pickle.load(handle)
+        except Exception as exc:  # noqa: BLE001 - tolerate a corrupt delta
+            logging.warning("Skipping unreadable history delta %s: %s", hist_file, exc)
+            continue
+        if isinstance(part, dict):
+            full.update(part)
+    return full
+
+
 class CheckpointManager:
     """Handles serialization and reconstruction of the sampler state to allow exact deterministic restarts."""
 
@@ -25,56 +71,59 @@ def save(
         iteration: int,
         space_model: Any,
         scaler: Any,
-        bin_state: Dict[str, Any],
-        history: Dict[str, Any],
-        sampler_state: Dict[str, Any] | None = None,
+        bin_state: dict[str, Any],
+        history: dict[str, Any],
+        sampler_state: dict[str, Any] | None = None,
     ) -> None:
         """Save a complete state snapshot."""
         iter_dir = self.checkpoint_dir / f"iter_{iteration}"
         iter_dir.mkdir(exist_ok=True)
-        (iter_dir / "format_version").write_text(str(CHECKPOINT_FORMAT_VERSION))
 
         # 1. Save Space Model (TVAE or TICA)
         if space_model is not None:
-            with open(iter_dir / "space_model.pkl", "wb") as f:
-                pickle.dump(space_model, f)
+            _atomic_pickle(space_model, iter_dir / "space_model.pkl")
 
         if hasattr(space_model, "type"):
             if (
                 space_model.type in _TORCH_ENCODER_MODES
                 and getattr(space_model, "fitted", None) is not None
             ):
-                torch.save(space_model.fitted.state_dict(), iter_dir / "model.pt")
+                tmp = iter_dir / "model.pt.tmp"
+                torch.save(space_model.fitted.state_dict(), tmp)
+                os.replace(tmp, iter_dir / "model.pt")
             elif (
                 space_model.type == "tica"
                 and getattr(space_model, "model", None) is not None
             ):
-                with open(iter_dir / "model.pkl", "wb") as f:
-                    pickle.dump(space_model.model, f)
+                _atomic_pickle(space_model.model, iter_dir / "model.pkl")
 
         # 2. Save Feature Scaler
-        with open(iter_dir / "scaler.pkl", "wb") as f:
-            pickle.dump(scaler, f)
-            
+        _atomic_pickle(scaler, iter_dir / "scaler.pkl")
+
         # 3. Save Bins & Spawn History
-        with open(iter_dir / "bin_state.pkl", "wb") as f:
-            pickle.dump(bin_state, f)
-            
-        # Delta Checkpointing: Only save history since the last checkpoint
+        _atomic_pickle(bin_state, iter_dir / "bin_state.pkl")
+
+        # Delta checkpointing: each file stores only the history since the last
+        # checkpoint. load()/reconstruct_history() merge the deltas back into the
+        # full history. Writes are atomic so a crash can't truncate a delta and
+        # break the chain.
         last_ckpt = self._get_latest_checkpoint_before(iteration)
         delta_history = {
             k: v for k, v in history.items() if k > last_ckpt and k <= iteration
         }
-        with open(iter_dir / "history.pkl", "wb") as f:
-            pickle.dump(delta_history, f)
-            
+        _atomic_pickle(delta_history, iter_dir / "history.pkl")
+
         if sampler_state is not None:
-            with open(iter_dir / "sampler_state.pkl", "wb") as f:
-                pickle.dump(sampler_state, f)
+            _atomic_pickle(sampler_state, iter_dir / "sampler_state.pkl")
+
+        # Write the format marker last: its presence signals a complete checkpoint.
+        tmp = iter_dir / "format_version.tmp"
+        tmp.write_text(str(CHECKPOINT_FORMAT_VERSION))
+        os.replace(tmp, iter_dir / "format_version")
 
     def load(
         self, iteration: int, space_model: Any = None
-    ) -> Tuple[Any, Any, Dict[str, Any], Dict[str, Any], Dict[str, Any]]:
+    ) -> tuple[Any, Any, dict[str, Any], dict[str, Any], dict[str, Any]]:
         """Restore the state exactly as it was at the specified iteration."""
         iter_dir = self.checkpoint_dir / f"iter_{iteration}"
         if not iter_dir.exists():
@@ -106,35 +155,12 @@ def load(
         # 2. Load scaler
         with open(iter_dir / "scaler.pkl", "rb") as f:
             scaler = pickle.load(f)
-            
-        # 3. Load bins & history
+
+        # 3. Load bins & reconstruct the full (delta-checkpointed) history.
         with open(iter_dir / "bin_state.pkl", "rb") as f:
             bin_state = pickle.load(f)
-            
-        with open(iter_dir / "history.pkl", "rb") as f:
-            history = pickle.load(f)
 
-        # Reconstruct full history for Delta Checkpointing
-        checkpoint_dirs = [
-            path
-            for path in self.checkpoint_dir.glob("iter_*")
-            if path.is_dir() and path.name.removeprefix("iter_").isdigit()
-        ]
-        previous_iters = sorted([
-            int(path.name.removeprefix("iter_"))
-            for path in checkpoint_dirs
-            if int(path.name.removeprefix("iter_")) < iteration
-        ], reverse=True)
-        
-        for prev_iter in previous_iters:
-            prev_hist_file = self.checkpoint_dir / f"iter_{prev_iter}" / "history.pkl"
-            if prev_hist_file.exists():
-                with open(prev_hist_file, "rb") as f:
-                    part_hist = pickle.load(f)
-                    if isinstance(part_hist, dict):
-                        for k, v in part_hist.items():
-                            if k not in history:
-                                history[k] = v
+        history = reconstruct_history(self.checkpoint_dir, iteration)
 
         state_path = iter_dir / "sampler_state.pkl"
         if state_path.exists():
diff --git a/autosampler/cli.py b/autosampler/cli.py
index d00aaff..7178a74 100644
--- a/autosampler/cli.py
+++ b/autosampler/cli.py
@@ -8,8 +8,9 @@
 import os
 import sys
 import tempfile
+from collections.abc import Sequence
 from pathlib import Path
-from typing import Any, Sequence
+from typing import Any
 
 os.environ.setdefault(
     "MPLCONFIGDIR", os.path.join(tempfile.gettempdir(), "autosampler-matplotlib")
diff --git a/autosampler/config.py b/autosampler/config.py
index 6c6d2a4..ef62a78 100644
--- a/autosampler/config.py
+++ b/autosampler/config.py
@@ -1,4 +1,4 @@
-from typing import Any, Dict, List, Optional, Union
+from typing import Any
 
 from pydantic import BaseModel, Field, field_validator, model_validator
 
@@ -7,10 +7,10 @@ class SystemConfig(BaseModel):
     conf_file: str
     top_file: str
     topology: str = "amber"
-    system_file: Optional[str] = None
-    project_file: Optional[str] = None
-    initial_trajectory: Optional[str] = None
-    trajectory_topology_file: Optional[str] = None
+    system_file: str | None = None
+    project_file: str | None = None
+    initial_trajectory: str | None = None
+    trajectory_topology_file: str | None = None
     feature_selection: str = "protein and not (type H)"
 
 
@@ -23,25 +23,25 @@ class EngineConfig(BaseModel):
     temperature: float = 300.0
     pressure: float = 1.0
     dt: float = 0.002
-    gpu_ids: Optional[List[int]] = None
-    gromacs_include_dir: Optional[str] = None
+    gpu_ids: list[int] | None = None
+    gromacs_include_dir: str | None = None
     gromacs_executable: str = "gmx"
-    gromacs_mdrun_nb: Optional[str] = None
-    gromacs_mdrun_pme: Optional[str] = None
-    gromacs_mdrun_update: Optional[str] = None
-    gromacs_mdrun_bonded: Optional[str] = None
-    gromacs_mdrun_pin: Optional[str] = None
+    gromacs_mdrun_nb: str | None = None
+    gromacs_mdrun_pme: str | None = None
+    gromacs_mdrun_update: str | None = None
+    gromacs_mdrun_bonded: str | None = None
+    gromacs_mdrun_pin: str | None = None
     gromacs_mdrun_ntmpi: int = 1
-    gromacs_mdrun_ntomp: Optional[int] = None
-    gromacs_mdrun_extra_args: List[str] = []
+    gromacs_mdrun_ntomp: int | None = None
+    gromacs_mdrun_extra_args: list[str] = []
     amber_executable: str = "pmemd"
-    amber_input_file: Optional[str] = None
-    amber_extra_args: List[str] = []
+    amber_input_file: str | None = None
+    amber_extra_args: list[str] = []
     amber_trajectory_format: str = "auto"
 
     @field_validator("gpu_ids")
     @classmethod
-    def validate_gpu_ids(cls, value: Optional[List[int]]) -> Optional[List[int]]:
+    def validate_gpu_ids(cls, value: list[int] | None) -> list[int] | None:
         if value is None:
             return None
         if not value:
@@ -67,12 +67,12 @@ class SpawningConfig(BaseModel):
     spawn_scheme: str = "density"
     spawn_type: str = "hard"
     search_mode: str = "explore"
-    n_bins: Optional[List[int]] = None
+    n_bins: list[int] | None = None
     walker: int = 10
     step: int = 10000
     stride: int = 100
     max_workers: int = 4
-    target: Optional[List[float]] = None
+    target: list[float] | None = None
     recent_density_window: int = 5
     voronoi_clusters: int = 150
     voronoi_periodic: bool = False
@@ -89,12 +89,12 @@ class AdaptiveModelConfig(BaseModel):
     lagtime: int = 5
     latent_dim: int = 2
     epochs: int = 50
-    batch_size: Union[int, str] = "auto"
+    batch_size: int | str = "auto"
     learning_rate: float = 0.0005
-    encoder_hidden_dims: List[int] = [256, 128]
-    decoder_hidden_dims: List[int] = [128, 256]
+    encoder_hidden_dims: list[int] = [256, 128]
+    decoder_hidden_dims: list[int] = [128, 256]
     dropout_rate: float = 0.1
-    deep_tica_hidden_dims: List[int] = [256, 128]
+    deep_tica_hidden_dims: list[int] = [256, 128]
     # SPIB (State Predictive Information Bottleneck) hyperparameters.
     spib_n_states: int = 10
     spib_beta: float = 1e-3
@@ -108,7 +108,7 @@ def validate_positive_int(cls, value: int) -> int:
 
     @field_validator("batch_size")
     @classmethod
-    def validate_batch_size(cls, value: Union[int, str]) -> Union[int, str]:
+    def validate_batch_size(cls, value: int | str) -> int | str:
         if isinstance(value, str):
             if value != "auto":
                 raise ValueError("batch_size must be 'auto' or a positive integer")
@@ -133,7 +133,7 @@ def validate_dropout_rate(cls, value: float) -> float:
 
     @field_validator("encoder_hidden_dims", "decoder_hidden_dims", "deep_tica_hidden_dims")
     @classmethod
-    def validate_hidden_dims(cls, value: List[int]) -> List[int]:
+    def validate_hidden_dims(cls, value: list[int]) -> list[int]:
         if not value:
             raise ValueError("hidden dimension lists must not be empty")
         if any(dim <= 0 for dim in value):
@@ -156,13 +156,13 @@ class MSMConfig(BaseModel):
     min_frames: int = 1000
     lagtime: int = 10
     # Optional lag-time ladder for an implied-timescale sweep (diagnostics).
-    lagtimes: Optional[List[int]] = None
+    lagtimes: list[int] | None = None
     n_microstates: int = 100
     cluster_method: str = "kmeans"  # "kmeans" | "regspace"
     estimator: str = "mle"  # "mle" | "bayesian"
     n_bayesian_samples: int = 50
     n_timescales: int = 3
-    n_metastable: Optional[int] = None
+    n_metastable: int | None = None
     # Keep microstate IDs comparable across iterations (seeds k-means from the
     # previous centres) — needed for transition-matrix convergence / spawning.
     stable_clustering: bool = False
@@ -171,7 +171,7 @@ class MSMConfig(BaseModel):
     spawn_leverage: int = 1  # slow eigenvectors used for the leverage factor
     spawn_uncertainty: bool = True  # include the outflow-uncertainty factor
     # Convergence: list of {name, params} criteria combined with all/any.
-    convergence_criteria: List[Dict[str, Any]] = Field(
+    convergence_criteria: list[dict[str, Any]] = Field(
         default_factory=lambda: [
             {"name": "implied_timescales", "params": {"tol": 0.1, "n_timescales": 2}},
             {"name": "vamp2", "params": {"tol": 0.05}},
@@ -225,12 +225,12 @@ class FeatureSelectionConfig(BaseModel):
     method: str = "greedy_vamp"  # "greedy_vamp" | "all"
     lagtime: int = 10
     cadence: int = 5  # re-select every N iterations (adaptive update)
-    max_features: Optional[int] = None  # cap on selected columns/groups
-    dim: Optional[int] = None  # singular values retained when scoring
+    max_features: int | None = None  # cap on selected columns/groups
+    dim: int | None = None  # singular values retained when scoring
     min_gain: float = 1e-4  # minimum VAMP-2 gain to add a feature group
     # Optional: rank these feature *types* by VAMP-2 and use the best one.
     # Empty -> always use the top-level `adaptive_feature_type`.
-    candidate_feature_types: List[str] = []
+    candidate_feature_types: list[str] = []
 
     @field_validator("method")
     @classmethod
@@ -241,7 +241,7 @@ def _method(cls, value: str) -> str:
 
     @field_validator("candidate_feature_types")
     @classmethod
-    def _candidate_types(cls, value: List[str]) -> List[str]:
+    def _candidate_types(cls, value: list[str]) -> list[str]:
         valid = {"distances", "fitted_coords", "phi_psi"}
         bad = [v for v in value if v not in valid]
         if bad:
@@ -269,19 +269,23 @@ class ExecutionConfig(BaseModel):
     """
 
     backend: str = "local"  # "local" | "slurm" | "pbs"
+    # Local backend: kill a walker (and abort the rest of the batch) if it runs
+    # longer than this many seconds. None disables the timeout (default). Guards
+    # against a hung in-process OpenMM walker stalling the campaign forever.
+    walker_timeout: float | None = None
     # Scheduler resource requests (per array task = one walker).
-    partition: Optional[str] = None  # SLURM partition / PBS queue
-    account: Optional[str] = None
+    partition: str | None = None  # SLURM partition / PBS queue
+    account: str | None = None
     walltime: str = "01:00:00"
     cpus_per_task: int = 1
     gpus_per_task: int = 0
-    memory: Optional[str] = None  # e.g. "8G"
+    memory: str | None = None  # e.g. "8G"
     # Robustness / polling.
     max_retries: int = 1  # resubmit failed walkers up to this many times
     poll_interval: float = 30.0  # seconds between scheduler polls
     submit_timeout: float = 60.0  # seconds for a submit/poll command
-    module_loads: List[str] = []  # `module load ...` lines for job scripts
-    extra_directives: List[str] = []  # raw #SBATCH / #PBS lines
+    module_loads: list[str] = []  # `module load ...` lines for job scripts
+    extra_directives: list[str] = []  # raw #SBATCH / #PBS lines
     job_name: str = "autosampler"
 
     @field_validator("backend")
@@ -339,9 +343,9 @@ class AutoSamplerConfig(BaseModel):
         default_factory=FeatureSelectionConfig
     )
     space_mode: str = "fixed"
-    n_bins: List[int] = [30, 30]
-    min_values: Optional[List[float]] = None
-    max_values: Optional[List[float]] = None
+    n_bins: list[int] = [30, 30]
+    min_values: list[float] | None = None
+    max_values: list[float] | None = None
     outdir: str = "runs/sampler_output"
     random_seed: int = 42
     checkpoint_freq: int = 1
@@ -352,7 +356,7 @@ class AutoSamplerConfig(BaseModel):
     retrain_policy: str = "fixed"
     vamp_retrain_tol: float = 0.1
     retrain_min_interval: int = 1
-    retrain_max_interval: Optional[int] = None
+    retrain_max_interval: int | None = None
     aggregate_memory: bool = True
     max_adaptive_memory_frames: int = 50000
     adaptive_feature_type: str = "distances"
@@ -377,7 +381,7 @@ def validate_space_mode(cls, value: str) -> str:
 
     @model_validator(mode="before")
     @classmethod
-    def promote_spawning_n_bins(cls, values: Dict[str, Any]) -> Dict[str, Any]:
+    def promote_spawning_n_bins(cls, values: dict[str, Any]) -> dict[str, Any]:
         values = dict(values)
         if "n_bins" not in values:
             spawning = values.get("spawning")
diff --git a/autosampler/core.py b/autosampler/core.py
index 5b1f348..9039ef5 100644
--- a/autosampler/core.py
+++ b/autosampler/core.py
@@ -1,10 +1,10 @@
-import logging
-import warnings
-from pathlib import Path
 import importlib.util
 import json
+import logging
 import shutil
 import sys
+import warnings
+from pathlib import Path
 
 # Suppress common non-critical warnings from dependencies
 warnings.filterwarnings(
@@ -12,7 +12,7 @@
 )
 warnings.filterwarnings("ignore", message="Reload offsets from trajectory")
 warnings.filterwarnings("ignore", message=".*Reader has no dt information.*")
-from typing import Any, Dict, List
+from typing import Any
 
 import numpy as np
 from pydantic import ValidationError
@@ -33,7 +33,7 @@
 class AutoSamplerCore:
     """Main orchestrator for the AutoSampler framework."""
 
-    def __init__(self, config_dict: Dict[str, Any]):
+    def __init__(self, config_dict: dict[str, Any]):
         try:
             self.config = AutoSamplerConfig(**config_dict)
         except ValidationError as e:
@@ -206,6 +206,7 @@ def validate_preflight(self) -> None:
     def _adaptive_model_kwargs(self) -> dict:
         kwargs = self.config.adaptive_model.model_dump()
         kwargs["space_mode"] = self.config.space_mode
+        kwargs["seed"] = self.config.random_seed  # reproducible CV training
         return kwargs
 
     def restore_checkpoint(self, iteration: int):
@@ -240,7 +241,7 @@ def restore_checkpoint(self, iteration: int):
     def latest_checkpoint_iteration(self) -> int:
         return self.checkpoint_manager.latest_iteration()
 
-    def resume_walkers(self) -> List[Any]:
+    def resume_walkers(self) -> list[Any]:
         """Rebuild next walkers from the latest restored checkpoint history."""
         if not self.history:
             return [self.engine.positions for _ in range(self.config.spawning.walker)]
@@ -278,38 +279,40 @@ def resume_walkers(self) -> List[Any]:
         )
 
 
-    def generate_initial_walkers(self) -> List[Any]:
+    def generate_initial_walkers(self) -> list[Any]:
         if not self.config.system.initial_trajectory:
             return [self.engine.positions for _ in range(self.config.spawning.walker)]
-            
+
         import logging
-        import numpy as np
         import random
         from pathlib import Path
+
         import MDAnalysis as mda
+        import numpy as np
+
         from autosampler.spaces import FeatureExtractor
-        
+
         traj_path = str(Path(self.config.system.initial_trajectory).resolve())
         logging.info(f"Initializing walkers from trajectory: {traj_path}")
-        
+
         trajectory_topology = (
             self.config.system.trajectory_topology_file or self.config.system.top_file
         )
-        
+
         u = mda.Universe(trajectory_topology, traj_path)
         n_frames = len(u.trajectory)
         n_walkers = self.config.spawning.walker
-        
+
         if n_frames == 0:
             raise ValueError(f"Initial trajectory {traj_path} contains 0 frames.")
-            
+
         points = None
         if hasattr(self, "_extract_physical_cvs") and self.config.system.project_file:
             try:
                 points = self._extract_physical_cvs([traj_path])
             except Exception as e:
                 logging.warning(f"Failed to extract CVs from initial trajectory, falling back to random sampling: {e}")
-                
+
         if points is not None and len(points) > n_walkers:
             logging.info(f"Selecting {n_walkers} starting walkers from {n_frames} frames using spawning scheme...")
             try:
@@ -327,7 +330,7 @@ def generate_initial_walkers(self) -> List[Any]:
             else:
                 logging.info(f"Only {n_frames} frames available for {n_walkers} walkers; replicating randomly.")
                 spawn_indices = [random.choice(range(n_frames)) for _ in range(n_walkers)]
-                
+
         feature_extractor = FeatureExtractor(
             topology=trajectory_topology,
             selection=self.config.system.feature_selection,
@@ -335,7 +338,7 @@ def generate_initial_walkers(self) -> List[Any]:
         walkers = feature_extractor.extract_positions_by_indices(
             [traj_path], spawn_indices
         )
-        
+
         if points is not None:
             from autosampler.core import build_frame_records
             frames = build_frame_records(
@@ -345,9 +348,9 @@ def generate_initial_walkers(self) -> List[Any]:
                 walker_parents=["initial"],
                 expected_frames=n_frames,
             )
-            
+
             next_walker_parents = [frames[idx]["key"] for idx in spawn_indices]
-            
+
             self.history[-1] = {
                 "projection": points,
                 "spawning_scheme": "initial",
@@ -359,7 +362,7 @@ def generate_initial_walkers(self) -> List[Any]:
             }
             self.walker_parents = next_walker_parents
             logging.info(f"Injected {n_frames} frames from initial trajectory into permanent history (iteration -1).")
-            
+
         return walkers
 
     def _traj_suffix(self) -> str:
@@ -370,7 +373,7 @@ def _traj_suffix(self) -> str:
             )
         return "xtc"
 
-    def run_iteration(self, walkers: List[Any]):
+    def run_iteration(self, walkers: list[Any]):
         """Run a single adaptive sampling iteration."""
 
         # 1. Run production MD
@@ -690,8 +693,8 @@ def run_iteration(self, walkers: List[Any]):
             "convergence_reason": self.convergence_reason,
         }
 
-    def _sampling_trajectories(self, current_trajectories: List[str]) -> List[str]:
-        trajectories: List[str] = []
+    def _sampling_trajectories(self, current_trajectories: list[str]) -> list[str]:
+        trajectories: list[str] = []
         for iteration in sorted(self.history):
             entry = self.history[iteration]
             if not isinstance(entry, dict) or entry.get("projection") is None:
@@ -707,9 +710,9 @@ def _sampling_trajectories(self, current_trajectories: List[str]) -> List[str]:
         return trajectories
 
     def _sampling_frame_records(
-        self, current_frame_records: List[Dict[str, Any]]
-    ) -> List[Dict[str, Any]]:
-        records: List[Dict[str, Any]] = []
+        self, current_frame_records: list[dict[str, Any]]
+    ) -> list[dict[str, Any]]:
+        records: list[dict[str, Any]] = []
         for iteration in sorted(self.history):
             entry = self.history[iteration]
             if not isinstance(entry, dict) or entry.get("projection") is None:
@@ -734,7 +737,7 @@ def _sampling_frame_records(
 
     def _infer_iteration_trajectories(
         self, iteration: int, projection: Any
-    ) -> List[str]:
+    ) -> list[str]:
         projection = np.asarray(projection)
         frames_per_walker = self.config.spawning.step // self.config.spawning.stride
         if frames_per_walker <= 0:
@@ -800,7 +803,7 @@ def _restore_walker_parents_from_history(self) -> None:
         if isinstance(latest_entry, dict):
             self.walker_parents = list(latest_entry.get("next_walker_parents") or [])
 
-    def _checkpoint_state(self) -> Dict[str, Any]:
+    def _checkpoint_state(self) -> dict[str, Any]:
         return {
             "n_bins": list(self.config.n_bins),
             "voronoi_clusters": self.config.spawning.voronoi_clusters,
@@ -817,7 +820,7 @@ def _checkpoint_state(self) -> Dict[str, Any]:
             "retrain_controller": self.retrain_controller.state_dict(),
         }
 
-    def _restore_sampler_state(self, state: Dict[str, Any] | None) -> None:
+    def _restore_sampler_state(self, state: dict[str, Any] | None) -> None:
         state = state or {}
         if "n_bins" in state:
             self.config.n_bins = list(state["n_bins"])
@@ -841,7 +844,7 @@ def _restore_sampler_state(self, state: Dict[str, Any] | None) -> None:
             self.msm_monitor.load_state_dict(state["msm_monitor"])
 
     @staticmethod
-    def _trajectory_file_problems(trajectories: List[str]) -> list[str]:
+    def _trajectory_file_problems(trajectories: list[str]) -> list[str]:
         bad: list[str] = []
         for path in trajectories:
             p = Path(path)
@@ -852,7 +855,7 @@ def _trajectory_file_problems(trajectories: List[str]) -> list[str]:
         return bad
 
     @staticmethod
-    def _validate_trajectory_files(trajectories: List[str]) -> None:
+    def _validate_trajectory_files(trajectories: list[str]) -> None:
         """Ensure each expected trajectory exists and is non-empty before reading.
 
         A walker can report success yet leave a missing or truncated file (disk
@@ -868,7 +871,7 @@ def _validate_trajectory_files(trajectories: List[str]) -> None:
 
     @staticmethod
     def _validate_sampling_trajectories(
-        trajectories: List[str], context: str = "sampling"
+        trajectories: list[str], context: str = "sampling"
     ) -> None:
         """Ensure cumulative trajectories are still available before spawning.
 
@@ -920,7 +923,7 @@ def _fs_apply(self, features: np.ndarray) -> np.ndarray:
         return np.asarray(features)[:, self.feature_selection_indices]
 
     def _extract_feature_type(
-        self, feature_extractor: FeatureExtractor, trajectories: List[str], ftype: str
+        self, feature_extractor: FeatureExtractor, trajectories: list[str], ftype: str
     ) -> np.ndarray:
         if ftype == "phi_psi":
             return feature_extractor.extract_aib9_phi_psi(trajectories)
@@ -929,7 +932,7 @@ def _extract_feature_type(
         return feature_extractor.extract_pairwise_distances(trajectories)
 
     def _extract_adaptive_features(
-        self, feature_extractor: FeatureExtractor, trajectories: List[str]
+        self, feature_extractor: FeatureExtractor, trajectories: list[str]
     ) -> np.ndarray:
         """Extract input features, optionally ranking candidate feature *types* by VAMP-2."""
         fs_cfg = getattr(self.config, "feature_selection", None)
@@ -954,7 +957,7 @@ def _extract_adaptive_features(
         from autosampler.spaces.feature_selection import rank_candidates
 
         n_frames = self.config.spawning.step // self.config.spawning.stride
-        extracted: Dict[str, np.ndarray] = {}
+        extracted: dict[str, np.ndarray] = {}
         for ftype in candidates:
             try:
                 extracted[ftype] = self._extract_feature_type(
@@ -1046,7 +1049,7 @@ def _maybe_select_features(
             selection.score,
         )
 
-    def _extract_physical_cvs(self, trajectories: List[str]) -> np.ndarray:
+    def _extract_physical_cvs(self, trajectories: list[str]) -> np.ndarray:
         """Load the user project file and extract physical CVs for ``trajectories``.
 
         Centralises the ``project_file`` import + ``extract_cvs`` call that was
@@ -1066,7 +1069,7 @@ def _extract_physical_cvs(self, trajectories: List[str]) -> np.ndarray:
             conf_file=self.config.system.conf_file,
         )
 
-    def _collect_msm_trajectories(self) -> List[Any]:
+    def _collect_msm_trajectories(self) -> list[Any]:
         """Split cumulative history projections into continuous per-walker trajectories.
 
         Each short walker is one continuous trajectory; transition counts are
@@ -1074,7 +1077,7 @@ def _collect_msm_trajectories(self) -> List[Any]:
         dropped because they carry no transitions.
         """
         frames_per_walker = self.config.spawning.step // self.config.spawning.stride
-        trajs: List[Any] = []
+        trajs: list[Any] = []
         for iteration in sorted(self.history):
             entry = self.history[iteration]
             if not isinstance(entry, dict):
@@ -1245,7 +1248,7 @@ def _increase_sampling_resolution(self) -> bool:
         new_bins = [min(int(b * 1.15), max_bins) for b in old_bins]
         new_bins = [
             min(old + 1, max_bins) if new <= old and old < max_bins else new
-            for old, new in zip(old_bins, new_bins)
+            for old, new in zip(old_bins, new_bins, strict=False)
         ]
         if new_bins == old_bins:
             return False
@@ -1450,13 +1453,13 @@ def _append_iteration_log(
         iteration: int,
         runner_time: float,
         other_time: float,
-        success: List[bool],
+        success: list[bool],
         points: Any,
         occupied_bins: int | None,
         total_bins: int | None,
         cumulative_frames: int | None,
-        spawn_indices: List[int],
-        trajectories: List[str],
+        spawn_indices: list[int],
+        trajectories: list[str],
     ) -> None:
         points_array = np.asarray(points)
         frames_this_iteration = len(points_array)
diff --git a/autosampler/engines/__init__.py b/autosampler/engines/__init__.py
index 8875240..6aa4b80 100644
--- a/autosampler/engines/__init__.py
+++ b/autosampler/engines/__init__.py
@@ -1,8 +1,31 @@
-from .amber import AmberEngine
+"""MD engines.
+
+Engine backends are registered *lazily*: each engine module (and its heavy,
+optional dependency — OpenMM / GROMACS / Amber) is imported only when that engine
+is first requested via ``EngineFactory.get``. This keeps ``import autosampler``
+free of the MD backends, so the base ``pip install autosampler`` does not need
+OpenMM (install it with ``pip install 'autosampler[openmm]'`` or via conda).
+"""
+
 from .base import EngineFactory, MDEngine
-from .gromacs import GromacsEngine
-from .openmm import OpenMMEngine
 
-EngineFactory.register("openmm", OpenMMEngine)
-EngineFactory.register("amber", AmberEngine)
-EngineFactory.register("gromacs", GromacsEngine)
+EngineFactory.register_lazy("openmm", "autosampler.engines.openmm", "OpenMMEngine")
+EngineFactory.register_lazy("amber", "autosampler.engines.amber", "AmberEngine")
+EngineFactory.register_lazy("gromacs", "autosampler.engines.gromacs", "GromacsEngine")
+
+__all__ = ["EngineFactory", "MDEngine"]
+
+
+def __getattr__(name):
+    # Backwards-compatible lazy access to the engine classes.
+    _classes = {
+        "OpenMMEngine": ("autosampler.engines.openmm", "OpenMMEngine"),
+        "AmberEngine": ("autosampler.engines.amber", "AmberEngine"),
+        "GromacsEngine": ("autosampler.engines.gromacs", "GromacsEngine"),
+    }
+    if name in _classes:
+        import importlib
+
+        module_path, class_name = _classes[name]
+        return getattr(importlib.import_module(module_path), class_name)
+    raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
diff --git a/autosampler/engines/amber.py b/autosampler/engines/amber.py
index 0f67ac0..2c28d4f 100644
--- a/autosampler/engines/amber.py
+++ b/autosampler/engines/amber.py
@@ -22,7 +22,6 @@
 import subprocess
 from functools import lru_cache
 from pathlib import Path
-from typing import Optional
 
 import numpy as np
 
@@ -59,7 +58,7 @@ def amber_trajectory_suffix(
 
 
 @lru_cache(maxsize=1)
-def _find_libnvjitlink_dir() -> Optional[str]:
+def _find_libnvjitlink_dir() -> str | None:
     """Find the directory containing libnvJitLink.so.12 dynamically."""
     import sys
 
@@ -139,8 +138,8 @@ def __init__(
         dt: float = 0.002,
         npt: bool = False,
         amber_executable: str = "pmemd",
-        amber_input_file: Optional[str] = None,
-        amber_extra_args: Optional[list[str]] = None,
+        amber_input_file: str | None = None,
+        amber_extra_args: list[str] | None = None,
         amber_trajectory_format: str = "auto",
         **kwargs,  # absorb OpenMM-specific kwargs (precision, platform_name, …)
     ):
@@ -155,16 +154,16 @@ def __init__(
         self._resolved_trajectory_format()
 
         # Set after prepare()
-        self.topology_file: Optional[Path] = None
-        self.start_coords_file: Optional[Path] = None
-        self.positions: Optional[str] = None  # str path for first-iteration walkers
+        self.topology_file: Path | None = None
+        self.start_coords_file: Path | None = None
+        self.positions: str | None = None  # str path for first-iteration walkers
 
     # ------------------------------------------------------------------
     # MDEngine interface
     # ------------------------------------------------------------------
 
     def prepare(
-        self, conf: Path, top: Path, system_file: Optional[Path] = None
+        self, conf: Path, top: Path, system_file: Path | None = None
     ) -> None:
         """Validate the Amber topology and coordinate files.
 
@@ -370,7 +369,7 @@ def _write_input(
         filepath: str,
         steps: int,
         stride: int,
-        trajectory_format: Optional[str] = None,
+        trajectory_format: str | None = None,
     ) -> None:
         """Write an Amber MD input (``.in``) file.
 
@@ -468,7 +467,7 @@ def _convert_nc_to_xtc(self, nc_file: str, xtc_file: str) -> None:
     def _write_rst7(
         filepath: str,
         positions_ang: np.ndarray,
-        box_ang: Optional[np.ndarray] = None,
+        box_ang: np.ndarray | None = None,
     ) -> None:
         """Write an Amber RST7 restart file (coordinates only, no velocities).
 
@@ -536,7 +535,7 @@ def _openmm_positions_to_angstrom(positions) -> np.ndarray:
         )
 
     @staticmethod
-    def _openmm_box_to_angstrom(box_vectors) -> Optional[np.ndarray]:
+    def _openmm_box_to_angstrom(box_vectors) -> np.ndarray | None:
         """Convert an OpenMM box-vectors tuple to ``[a, b, c, 90, 90, 90]`` in Angstroms.
 
         Only orthogonal boxes are currently supported; triclinic cells will
diff --git a/autosampler/engines/base.py b/autosampler/engines/base.py
index a074b0b..4c2d3d2 100644
--- a/autosampler/engines/base.py
+++ b/autosampler/engines/base.py
@@ -37,15 +37,42 @@ def run_production(self, run_index: int, start_coords: Path, steps: int,
 # Factory Registry
 class EngineFactory:
     _engines = {}
+    _lazy: dict = {}  # name -> (module_path, class_name)
 
     @classmethod
     def register(cls, name: str, engine_cls):
         """Register a new engine implementation."""
         cls._engines[name] = engine_cls
 
+    @classmethod
+    def register_lazy(cls, name: str, module_path: str, class_name: str):
+        """Register an engine that is imported only when first requested.
+
+        Keeps heavy optional backends (OpenMM, GROMACS, Amber) out of the import
+        path of ``import autosampler`` so the base install need not pull them in.
+        """
+        cls._lazy[name] = (module_path, class_name)
+
     @classmethod
     def get(cls, name: str, **kwargs) -> MDEngine:
-        """Instantiate an engine by name."""
+        """Instantiate an engine by name (importing its backend on first use)."""
+        if name not in cls._engines and name in cls._lazy:
+            import importlib
+
+            module_path, class_name = cls._lazy[name]
+            try:
+                module = importlib.import_module(module_path)
+            except ImportError as exc:
+                raise ImportError(
+                    f"MD engine {name!r} needs an optional dependency that is not "
+                    f"installed ({exc}). Install it, e.g. "
+                    f"`pip install 'autosampler[{name}]'` or via conda."
+                ) from exc
+            cls._engines[name] = getattr(module, class_name)
         if name not in cls._engines:
             raise ValueError(f"Unknown MD engine: {name}")
         return cls._engines[name](**kwargs)
+
+    @classmethod
+    def available(cls) -> list:
+        return sorted(set(cls._engines) | set(cls._lazy))
diff --git a/autosampler/engines/gromacs.py b/autosampler/engines/gromacs.py
index 41bc9b4..7658857 100644
--- a/autosampler/engines/gromacs.py
+++ b/autosampler/engines/gromacs.py
@@ -28,7 +28,6 @@
 import shutil
 import subprocess
 from pathlib import Path
-from typing import Optional
 
 import numpy as np
 
@@ -65,15 +64,15 @@ def __init__(
         dt: float = 0.002,
         npt: bool = False,
         gromacs_executable: str = "gmx",
-        gromacs_include_dir: Optional[str] = None,
-        gromacs_mdrun_nb: Optional[str] = None,
-        gromacs_mdrun_pme: Optional[str] = None,
-        gromacs_mdrun_update: Optional[str] = None,
-        gromacs_mdrun_bonded: Optional[str] = None,
-        gromacs_mdrun_pin: Optional[str] = None,
+        gromacs_include_dir: str | None = None,
+        gromacs_mdrun_nb: str | None = None,
+        gromacs_mdrun_pme: str | None = None,
+        gromacs_mdrun_update: str | None = None,
+        gromacs_mdrun_bonded: str | None = None,
+        gromacs_mdrun_pin: str | None = None,
         gromacs_mdrun_ntmpi: int = 1,
-        gromacs_mdrun_ntomp: Optional[int] = None,
-        gromacs_mdrun_extra_args: Optional[list[str]] = None,
+        gromacs_mdrun_ntomp: int | None = None,
+        gromacs_mdrun_extra_args: list[str] | None = None,
         **kwargs,  # absorb OpenMM / Amber specific kwargs
     ):
         self.temperature = temperature
@@ -92,17 +91,17 @@ def __init__(
         self.gromacs_mdrun_extra_args = list(gromacs_mdrun_extra_args or [])
 
         # Set after prepare()
-        self.topology_file: Optional[Path] = None
-        self.start_coords_file: Optional[Path] = None
-        self.mdp_template: Optional[Path] = None
-        self.positions: Optional[str] = None  # str path for first-iteration walkers
+        self.topology_file: Path | None = None
+        self.start_coords_file: Path | None = None
+        self.mdp_template: Path | None = None
+        self.positions: str | None = None  # str path for first-iteration walkers
 
     # ------------------------------------------------------------------
     # MDEngine interface
     # ------------------------------------------------------------------
 
     def prepare(
-        self, conf: Path, top: Path, system_file: Optional[Path] = None
+        self, conf: Path, top: Path, system_file: Path | None = None
     ) -> None:
         """Validate GROMACS input files and store run-time parameters.
 
@@ -448,7 +447,7 @@ def _write_gro(
         self,
         filepath: str,
         positions_ang: np.ndarray,
-        box_ang: Optional[np.ndarray] = None,
+        box_ang: np.ndarray | None = None,
     ) -> None:
         """Write a GROMACS GRO file with updated positions.
 
@@ -518,7 +517,7 @@ def _openmm_positions_to_angstrom(positions) -> np.ndarray:
         )
 
     @staticmethod
-    def _openmm_box_to_angstrom(box_vectors) -> Optional[np.ndarray]:
+    def _openmm_box_to_angstrom(box_vectors) -> np.ndarray | None:
         """Convert an OpenMM box-vectors tuple to ``[a, b, c, 90, 90, 90]`` in Angstroms.
 
         Only orthogonal boxes are currently supported.
diff --git a/autosampler/engines/openmm.py b/autosampler/engines/openmm.py
index a34dcbb..1e5917d 100644
--- a/autosampler/engines/openmm.py
+++ b/autosampler/engines/openmm.py
@@ -2,7 +2,6 @@
 import inspect
 import os
 from pathlib import Path
-from typing import Optional
 
 from openmm import *
 from openmm.app import *
@@ -23,7 +22,7 @@ def __init__(
         precision: str = "mixed",
         npt: bool = False,
         equilibrate: bool = False,
-        gromacs_include_dir: Optional[str] = None,
+        gromacs_include_dir: str | None = None,
         **kwargs,
     ):
         self.temperature_val = temperature
@@ -50,20 +49,13 @@ def _get_platform(cls, platform_name: str):
         try:
             return Platform.getPlatformByName(platform_name)
         except Exception as exc:
-            available = ", ".join(cls._available_platforms()) or "none"
-            hint = (
-                "Install the CUDA-enabled OpenMM package, for example "
-                "`conda install -c conda-forge openmm cuda-version=12` or "
-                "`python -m pip install 'openmm[cuda12]'`, after confirming "
-                "that the NVIDIA driver is installed. For CPU-only runs, set "
-                "`engine.platform_name: CPU` in the AutoSampler config."
-            )
+            ", ".join(cls._available_platforms()) or "none"
             import logging
             logging.warning(f"OpenMM platform {platform_name} validation failed (likely because you are on a login node). Assuming compute nodes will have it. Error: {exc}")
             return None
 
     def prepare(
-        self, conf: Path, top: Path, system_file: Optional[Path] = None
+        self, conf: Path, top: Path, system_file: Path | None = None
     ) -> None:
         """Prepare the MD environment, e.g., setup system, topology, forces."""
         gro_file = str(conf)
diff --git a/autosampler/execution/__init__.py b/autosampler/execution/__init__.py
index 8efa5de..81ed50d 100644
--- a/autosampler/execution/__init__.py
+++ b/autosampler/execution/__init__.py
@@ -39,7 +39,10 @@ def make_backend(execution_config, *, gpu_ids=None, max_workers: int = 8):
     backend = getattr(execution_config, "backend", "local")
     if backend == "local":
         return ExecutionBackendFactory.get(
-            "local", gpu_ids=gpu_ids, max_workers=max_workers
+            "local",
+            gpu_ids=gpu_ids,
+            max_workers=max_workers,
+            walker_timeout=getattr(execution_config, "walker_timeout", None),
         )
 
     cfg = execution_config.model_dump()
diff --git a/autosampler/execution/local.py b/autosampler/execution/local.py
index e39b43e..4759fa7 100644
--- a/autosampler/execution/local.py
+++ b/autosampler/execution/local.py
@@ -72,18 +72,47 @@ def _execution_slots(
 
 def _run_one(task: WalkerTask, device_index: int) -> bool:
     task.device_index = device_index
-    return run_walker_task(task)
+    # A single walker's failure (CUDA error, NaN blow-up, missing file, …) must
+    # not abort the whole iteration — mirror the scheduler path (run_task.py),
+    # which reports failures as success=False rather than raising.
+    try:
+        return run_walker_task(task)
+    except Exception:  # noqa: BLE001 - report failure, keep the batch alive
+        import logging
+        import traceback
+
+        logging.error(
+            "Walker %s failed; marking it unsuccessful and continuing:\n%s",
+            getattr(task, "index", "?"),
+            traceback.format_exc(),
+        )
+        return False
+
+
+def _terminate_workers(executor) -> None:
+    """Best-effort kill of a ProcessPoolExecutor's worker processes."""
+    for proc in list(getattr(executor, "_processes", {}).values()):
+        try:
+            proc.terminate()
+        except Exception:  # noqa: BLE001 - already gone / not terminable
+            pass
 
 
 class LocalProcessBackend(ExecutionBackend):
     def __init__(
-        self, gpu_ids: list[int] | None = None, max_workers: int = 8, **_
+        self,
+        gpu_ids: list[int] | None = None,
+        max_workers: int = 8,
+        walker_timeout: float | None = None,
+        **_,
     ):
         self.gpu_ids = gpu_ids
         self.max_workers = max_workers
+        self.walker_timeout = walker_timeout
 
     def execute(self, tasks: list[WalkerTask]) -> list[bool]:
         import multiprocessing as mp
+        import time
 
         if not tasks:
             return []
@@ -97,6 +126,7 @@ def execute(self, tasks: list[WalkerTask]) -> list[bool]:
         ctx = mp.get_context("spawn")
         results = [False] * len(tasks)
         task_iter = iter(tasks)
+        timeout = self.walker_timeout
 
         def submit(executor, device_index: int):
             try:
@@ -104,25 +134,56 @@ def submit(executor, device_index: int):
             except StopIteration:
                 return None
             future = executor.submit(_run_one, task, device_index)
-            return future, task.index, device_index
+            # value = (index, device, start_time)
+            return future, (task.index, device_index, time.monotonic())
 
         with ProcessPoolExecutor(max_workers=len(slots), mp_context=ctx) as executor:
             active: dict = {}
             for device_index in slots:
                 submitted = submit(executor, device_index)
                 if submitted is not None:
-                    future, idx, dev = submitted
-                    active[future] = (idx, dev)
+                    future, meta = submitted
+                    active[future] = meta
 
             while active:
-                done, _ = wait(active, return_when=FIRST_COMPLETED)
+                # Poll at the timeout cadence so overdue walkers are detected even
+                # when nothing completes; `wait` itself returns no timed-out futures.
+                poll = None if timeout is None else max(min(timeout, 30.0), 0.05)
+                done, _ = wait(active, timeout=poll, return_when=FIRST_COMPLETED)
                 for future in done:
-                    idx, freed_device = active.pop(future)
-                    results[idx] = future.result()
+                    idx, freed_device, _ = active.pop(future)
+                    try:
+                        results[idx] = future.result()
+                    except Exception:  # noqa: BLE001 - e.g. a killed worker process
+                        import logging
+
+                        logging.error(
+                            "Walker %s did not return a result (worker died); "
+                            "marking it unsuccessful.",
+                            idx,
+                        )
+                        results[idx] = False
                     submitted = submit(executor, freed_device)
                     if submitted is not None:
-                        next_future, next_idx, dev = submitted
-                        active[next_future] = (next_idx, dev)
+                        next_future, meta = submitted
+                        active[next_future] = meta
+
+                if timeout is not None and active:
+                    now = time.monotonic()
+                    overdue = [m for m in active.values() if now - m[2] > timeout]
+                    if overdue:
+                        import logging
+
+                        for idx, _dev, _start in active.values():
+                            results[idx] = False
+                        logging.error(
+                            "Walker(s) %s exceeded walker_timeout=%ss; terminating "
+                            "the batch and marking remaining walkers unsuccessful.",
+                            sorted(m[0] for m in overdue),
+                            timeout,
+                        )
+                        _terminate_workers(executor)
+                        active.clear()
 
         return results
 
diff --git a/autosampler/log_cli.py b/autosampler/log_cli.py
index a20a703..431d436 100644
--- a/autosampler/log_cli.py
+++ b/autosampler/log_cli.py
@@ -3,8 +3,8 @@
 from __future__ import annotations
 
 import argparse
+from collections.abc import Sequence
 from pathlib import Path
-from typing import Sequence
 
 from autosampler.cli import load_config
 from autosampler.logs import write_exploration_log
diff --git a/autosampler/path_cli.py b/autosampler/path_cli.py
index e0aad25..b0ccf07 100644
--- a/autosampler/path_cli.py
+++ b/autosampler/path_cli.py
@@ -7,9 +7,10 @@
 import json
 import re
 import sys
+from collections.abc import Sequence
 from dataclasses import dataclass
 from pathlib import Path
-from typing import Any, Sequence
+from typing import Any
 
 import numpy as np
 
diff --git a/autosampler/paths.py b/autosampler/paths.py
index dc35896..1a67b78 100644
--- a/autosampler/paths.py
+++ b/autosampler/paths.py
@@ -3,14 +3,14 @@
 from __future__ import annotations
 
 import json
-import pickle
 import re
 import shutil
 import subprocess
 import tempfile
+from collections.abc import Iterable
 from dataclasses import dataclass
 from pathlib import Path
-from typing import Any, Iterable
+from typing import Any
 
 import numpy as np
 
@@ -40,7 +40,7 @@ def to_dict(self) -> dict[str, Any]:
         }
 
     @classmethod
-    def from_dict(cls, data: dict[str, Any]) -> "FrameRef":
+    def from_dict(cls, data: dict[str, Any]) -> FrameRef:
         return cls(
             iteration=int(data["iteration"]),
             walker=int(data["walker"]),
@@ -97,7 +97,7 @@ def build_frame_records(
 
     records: list[dict[str, Any]] = []
     point_offset = 0
-    for walker, (trajectory, count) in enumerate(zip(trajectories, counts)):
+    for walker, (trajectory, count) in enumerate(zip(trajectories, counts, strict=False)):
         for frame in range(count):
             parent = (
                 walker_parents[walker]
@@ -126,26 +126,29 @@ def map_global_frame(records: list[dict[str, Any]], index: int) -> dict[str, Any
 
 
 def load_history(run_dir: Path, checkpoint: int | None = None) -> dict[int, Any]:
+    from autosampler.checkpoints.manager import reconstruct_history
+
     checkpoint_root = run_dir / "checkpoints"
     if checkpoint is None:
-        checkpoint_dirs = [
-            path
+        checkpoint_iters = [
+            int(path.name.removeprefix("iter_"))
             for path in checkpoint_root.glob("iter_*")
             if path.is_dir() and path.name.removeprefix("iter_").isdigit()
         ]
-        if not checkpoint_dirs:
+        if not checkpoint_iters:
             raise FileNotFoundError(f"No checkpoints found under {checkpoint_root}")
-        checkpoint_dir = max(
-            checkpoint_dirs, key=lambda path: int(path.name.removeprefix("iter_"))
-        )
+        target = max(checkpoint_iters)
     else:
-        checkpoint_dir = checkpoint_root / f"iter_{checkpoint}"
+        target = checkpoint
+        if not (checkpoint_root / f"iter_{target}").exists():
+            raise FileNotFoundError(
+                f"History file not found: {checkpoint_root / f'iter_{target}'}"
+            )
 
-    history_path = checkpoint_dir / "history.pkl"
-    if not history_path.exists():
-        raise FileNotFoundError(f"History file not found: {history_path}")
-    with history_path.open("rb") as handle:
-        return pickle.load(handle)
+    # History is delta-checkpointed: each iter_*/history.pkl holds only the
+    # entries since the previous checkpoint. Merge them back into the full
+    # history (otherwise the lineage/path tools see only the last window).
+    return reconstruct_history(checkpoint_root, target)
 
 
 def history_records(history: dict[int, Any]) -> list[FrameRef]:
diff --git a/autosampler/spaces/features.py b/autosampler/spaces/features.py
index ece9bb5..6670f74 100644
--- a/autosampler/spaces/features.py
+++ b/autosampler/spaces/features.py
@@ -1,4 +1,3 @@
-from typing import Optional
 
 import MDAnalysis as mda  # type: ignore
 import numpy as np
@@ -22,7 +21,7 @@ def _load_universe(topology: str, trajectories: list[str] | str, **kwargs) -> md
 class FeatureExtractor:
     """Extracts features from MD trajectories for dimensionality reduction."""
 
-    def __init__(self, topology: str, selection: Optional[str] = None):
+    def __init__(self, topology: str, selection: str | None = None):
         self.topology = topology
         self.selection = selection if selection else "protein and not (type H)"
 
@@ -32,14 +31,14 @@ def extract_pairwise_distances(self, trajectories: list[str]) -> np.ndarray:
         try:
             u = _load_universe(self.topology, trajectories)
         except Exception as e:
-            raise ValueError(f"Failed to load universe with top={self.topology}, trajs={trajectories}: {e}")
+            raise ValueError(f"Failed to load universe with top={self.topology}, trajs={trajectories}: {e}") from e
 
         ag = u.select_atoms(self.selection)
         num_pairs = ag.n_atoms * (ag.n_atoms - 1) // 2
 
         try:
             dist_list = np.zeros((u.trajectory.n_frames, num_pairs), dtype=np.float32)
-            for j, ts in enumerate(u.trajectory):
+            for j, _ts in enumerate(u.trajectory):
                 r = distance_array(ag, ag, box=u.dimensions, backend="OpenMP")
                 r = r[np.triu_indices(r.shape[0], k=1)]
                 dist_list[j] = r
@@ -57,14 +56,14 @@ def extract_fitted_coords(self, trajectories: list[str]) -> np.ndarray:
             ref = mda.Universe(self.topology)
             u = _load_universe(self.topology, trajectories)
         except Exception as e:
-            raise ValueError(f"Failed to load universe for fitted coords extraction: {e}")
+            raise ValueError(f"Failed to load universe for fitted coords extraction: {e}") from e
 
         ag = u.select_atoms(self.selection)
         num_features = ag.n_atoms * 3
 
         try:
             coord_list = np.zeros((u.trajectory.n_frames, num_features), dtype=np.float32)
-            for j, ts in enumerate(u.trajectory):
+            for j, _ts in enumerate(u.trajectory):
                 # Align the current frame in memory to the reference structure
                 align.alignto(u, ref, select=self.selection)
                 # Store the flattened coordinates of the selection
@@ -86,7 +85,7 @@ def extract_aib9_phi_psi(self, trajectories: list[str]) -> np.ndarray:
         try:
             u = _load_universe(self.topology, trajectories)
         except Exception as e:
-            raise ValueError(f"Failed to load universe for AIB9 phi/psi extraction: {e}")
+            raise ValueError(f"Failed to load universe for AIB9 phi/psi extraction: {e}") from e
 
         residues = list(u.select_atoms("resname AIB").residues)
         if len(residues) != 9:
@@ -116,7 +115,7 @@ def atom(residue, name: str):
             values = np.zeros((u.trajectory.n_frames, 18), dtype=np.float32)
             for frame_index, ts in enumerate(u.trajectory):
                 row = []
-                for phi_group, psi_group in zip(phi_atoms, psi_atoms):
+                for phi_group, psi_group in zip(phi_atoms, psi_atoms, strict=False):
                     phi = calc_dihedrals(
                         phi_group[0].positions,
                         phi_group[1].positions,
@@ -146,16 +145,16 @@ def extract_rg_rmsd(self, trajectories: list[str], reference_pdb: str) -> np.nda
             u = _load_universe(self.topology, trajectories)
             ref = _load_universe(self.topology, reference_pdb)
         except Exception as e:
-            raise ValueError(f"Failed to load universe for Rg/RMSD: {e}")
+            raise ValueError(f"Failed to load universe for Rg/RMSD: {e}") from e
 
         protein = u.select_atoms("protein")
 
         try:
             # Calculate Rg
             rg_values = np.zeros(u.trajectory.n_frames, dtype=np.float32)
-            for j, ts in enumerate(u.trajectory):
+            for j, _ts in enumerate(u.trajectory):
                 rg_values[j] = protein.radius_of_gyration()
-    
+
             # Calculate RMSD
             # Default to CA backbone for stable RMSD
             R = rms.RMSD(u, ref, select="protein", groupselections=["name CA"])
@@ -171,9 +170,8 @@ def extract_rg_rmsd(self, trajectories: list[str], reference_pdb: str) -> np.nda
 
     def extract_positions_by_indices(self, trajectories: list[str], indices: list[int]) -> list:
         """Return full-system OpenMM positions for selected global frame indices."""
-        from openmm import Vec3  # type: ignore
+        from MDAnalysis.coordinates.XTC import XTCReader  # type: ignore
         from openmm.unit import nanometer  # type: ignore
-        from MDAnalysis.coordinates.XTC import XTCReader # type: ignore
 
         # Map global indices to their original order to return them correctly
         sorted_requests = sorted(list(enumerate(indices)), key=lambda x: x[1])
@@ -187,11 +185,11 @@ def extract_positions_by_indices(self, trajectories: list[str], indices: list[in
             for original_order, target_index in sorted_requests:
                 if target_index < 0:
                     raise IndexError(f"Spawn frame index {target_index} cannot be negative.")
-                
+
                 while True:
                     if traj_idx >= len(trajectories):
                         raise IndexError(f"Spawn frame index {target_index} is outside trajectory range.")
-                        
+
                     # Compute length extremely fast without parsing topology
                     if str(trajectories[traj_idx]).endswith(".xtc"):
                         with XTCReader(trajectories[traj_idx]) as reader:
diff --git a/autosampler/spaces/model.py b/autosampler/spaces/model.py
index 60a15d2..7936492 100644
--- a/autosampler/spaces/model.py
+++ b/autosampler/spaces/model.py
@@ -1,12 +1,13 @@
-import torch
-import numpy as np
 from typing import Any
 
+import numpy as np
+import torch
 from deeptime.decomposition.deep import TVAE
 from deeptime.util.data import TrajectoryDataset
 from torch.utils.data import DataLoader
+
 from .scalers import TrajectoryScaler
-from .tvae import TVAEBottleneckEncoder, TVAEBottleneckDecoder
+from .tvae import TVAEBottleneckDecoder, TVAEBottleneckEncoder
 
 
 class AdaptiveSpaceModel:
@@ -24,6 +25,7 @@ class AdaptiveSpaceModel:
         "deep_tica_hidden_dims": [256, 128],
         "spib_n_states": 10,
         "spib_beta": 1e-3,
+        "seed": 0,
     }
 
     def __init__(
@@ -40,9 +42,11 @@ def __init__(
         deep_tica_hidden_dims: list[int] | None = None,
         spib_n_states: int = 10,
         spib_beta: float = 1e-3,
+        seed: int = 0,
         **_: Any,
     ):
         self.type = space_mode
+        self.seed = int(seed)
         self.lagtime = int(lagtime)
         self.latent_dim = int(latent_dim)
         self.epochs = int(epochs)
@@ -98,6 +102,13 @@ def fit(self, features: np.ndarray, walker_length: int, n_walkers: int):
         self.ensure_config_defaults()
         input_size = features.shape[-1]
 
+        # Reseed the torch RNG from the configured seed before every fit so that a
+        # CV retrained at iteration N is reproducible regardless of the RNG draws
+        # made in between (network init, DataLoader shuffling, etc.).
+        torch.manual_seed(self.seed)
+        if torch.cuda.is_available():
+            torch.cuda.manual_seed_all(self.seed)
+
         # Fail fast with an actionable message if an optional backend is missing.
         from .registry import ensure_available, is_adaptive_space
 
@@ -169,8 +180,8 @@ def fit(self, features: np.ndarray, walker_length: int, n_walkers: int):
             ]
             import lightning as pl
             from mlcolvar.cvs import DeepTICA
+            from mlcolvar.data import DictDataset, DictModule
             from mlcolvar.utils.timelagged import create_timelagged_dataset
-            from mlcolvar.data import DictModule, DictDataset
 
             data_dict = {}
             for traj in traj_list:
@@ -253,6 +264,7 @@ def fit(self, features: np.ndarray, walker_length: int, n_walkers: int):
                 beta=self.spib_beta,
                 dropout=self.dropout_rate,
                 device=self.device,
+                seed=self.seed,
             )
 
         elif self.type == "deep-lda":
@@ -299,15 +311,15 @@ def project(self, features: np.ndarray) -> np.ndarray:
                 mean, _ = self.fitted(tensor)
                 projected = mean.detach().cpu().numpy()
         elif self.type == "deep-tica":
+            try:
+                device = next(self.model.parameters()).device
+            except (StopIteration, AttributeError):
+                device = torch.device("cpu")
             tensor = torch.as_tensor(
-                self._torch_features(scaled), dtype=torch.float32
+                self._torch_features(scaled), dtype=torch.float32, device=device
             )
             with torch.no_grad():
-                projected = (
-                    self.model(tensor)
-                    .detach()
-                    .numpy()
-                )
+                projected = self.model(tensor).detach().cpu().numpy()
         elif self.type == "pca":
             projected = self.model.transform(scaled)
         elif self.type == "tica":
diff --git a/autosampler/spaces/scalers.py b/autosampler/spaces/scalers.py
index 7915541..3e50507 100644
--- a/autosampler/spaces/scalers.py
+++ b/autosampler/spaces/scalers.py
@@ -1,9 +1,10 @@
-from sklearn.preprocessing import MinMaxScaler, StandardScaler, RobustScaler
 import numpy as np
+from sklearn.preprocessing import MinMaxScaler, RobustScaler, StandardScaler
+
 
 class TrajectoryScaler:
     """Scales pairwise feature vectors while preserving shape mappings."""
-    
+
     def __init__(self, scaler_type: str = "minmax"):
         if scaler_type == "minmax":
             self.scaler = MinMaxScaler()
diff --git a/autosampler/spaces/tvae.py b/autosampler/spaces/tvae.py
index 78f6ee7..acb9d2b 100644
--- a/autosampler/spaces/tvae.py
+++ b/autosampler/spaces/tvae.py
@@ -1,12 +1,13 @@
+
 import torch
 import torch.nn as nn
-from typing import List, Tuple
+
 
 class TVAEBottleneckEncoder(nn.Module):
     """Bottleneck TVAE encoder with configurable hidden dimensions and dropout."""
-    
-    def __init__(self, input_dim: int, latent_dim: int = 2, 
-                 hidden_dims: List[int] | None = None, dropout_rate: float = 0.1):
+
+    def __init__(self, input_dim: int, latent_dim: int = 2,
+                 hidden_dims: list[int] | None = None, dropout_rate: float = 0.1):
         super().__init__()
         hidden_dims = hidden_dims or [256, 128]
         layers = []
@@ -19,20 +20,20 @@ def __init__(self, input_dim: int, latent_dim: int = 2,
                 layers.append(nn.Dropout(dropout_rate))
             curr_dim = h_dim
         self.backbone = nn.Sequential(*layers)
-        
+
         # Latent outputs (mean and log-variance)
         self.mean_layer = nn.Linear(curr_dim, latent_dim)
         self.logvar_layer = nn.Linear(curr_dim, latent_dim)
 
-    def forward(self, x: torch.Tensor) -> Tuple[torch.Tensor, torch.Tensor]:
+    def forward(self, x: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
         h = self.backbone(x)
         return self.mean_layer(h), self.logvar_layer(h)
 
 class TVAEBottleneckDecoder(nn.Module):
     """Bottleneck TVAE decoder to reconstruct the original features from the latent space."""
-    
-    def __init__(self, latent_dim: int, output_dim: int, 
-                 hidden_dims: List[int] | None = None, dropout_rate: float = 0.1):
+
+    def __init__(self, latent_dim: int, output_dim: int,
+                 hidden_dims: list[int] | None = None, dropout_rate: float = 0.1):
         super().__init__()
         hidden_dims = hidden_dims or [128, 256]
         layers = []
@@ -44,7 +45,7 @@ def __init__(self, latent_dim: int, output_dim: int,
             if dropout_rate > 0:
                 layers.append(nn.Dropout(dropout_rate))
             curr_dim = h_dim
-            
+
         self.backbone = nn.Sequential(*layers)
         self.output_layer = nn.Linear(curr_dim, output_dim)
 
diff --git a/autosampler/spawners/__init__.py b/autosampler/spawners/__init__.py
index 7668c13..b4a89b3 100644
--- a/autosampler/spawners/__init__.py
+++ b/autosampler/spawners/__init__.py
@@ -1,6 +1,6 @@
 from .base import Spawner, SpawnerFactory
-from .fps import FPSSpawner
 from .density import DensitySpawner
+from .fps import FPSSpawner
 from .lof import LOFSpawner
 from .msm import MSMSpawner
 from .voronoi import VoronoiSpawner
diff --git a/autosampler/spawners/base.py b/autosampler/spawners/base.py
index 19da8de..342a6ce 100644
--- a/autosampler/spawners/base.py
+++ b/autosampler/spawners/base.py
@@ -1,20 +1,22 @@
 from abc import ABC, abstractmethod
+from typing import Any
+
 import numpy as np
-from typing import Any, Dict, List, Optional
+
 
 class Spawner(ABC):
     """Abstract Strategy interface for spawning outliers from explored space."""
-    
+
     @abstractmethod
-    def sample(self, points: np.ndarray, top_n: int, history: Optional[Dict[int, Any]] = None) -> List[int]:
+    def sample(self, points: np.ndarray, top_n: int, history: dict[int, Any] | None = None) -> list[int]:
         """Select top_n points from the set of explored points.
-        
+
         Args:
             points: numpy array of shape (n_points, n_features)
             top_n: number of points to select
             history: optional sampler history used by spawners that need
                 cumulative state.
-            
+
         Returns:
             List of indices of the selected points.
         """
diff --git a/autosampler/spawners/density.py b/autosampler/spawners/density.py
index 56c4bba..0802a48 100644
--- a/autosampler/spawners/density.py
+++ b/autosampler/spawners/density.py
@@ -6,6 +6,7 @@
 import numpy as np
 
 from autosampler.binning.spatial import BinTable, RegularBinner
+
 from .base import Spawner, SpawnerFactory
 
 
@@ -124,7 +125,11 @@ def _weighted_choice(rows: np.ndarray, weights: np.ndarray, top_n: int) -> np.nd
     if weights.sum() <= 0:
         weights = np.ones_like(weights, dtype=float)
     weights = weights / weights.sum()
-    replace = len(rows) < top_n
+    # replace must be True when fewer than top_n entries have nonzero probability,
+    # otherwise np.random.choice raises "fewer non-zero entries in p than size"
+    # (common in target mode, where target_closeness zeroes out most bins).
+    n_nonzero = int(np.count_nonzero(weights))
+    replace = len(rows) < top_n or n_nonzero < top_n
     return np.random.choice(rows, size=top_n, replace=replace, p=weights).astype(int)
 
 
diff --git a/autosampler/spawners/fps.py b/autosampler/spawners/fps.py
index 7382bac..b3eff54 100644
--- a/autosampler/spawners/fps.py
+++ b/autosampler/spawners/fps.py
@@ -1,4 +1,3 @@
-from typing import Optional
 
 import numpy as np
 
@@ -17,7 +16,7 @@ class FPSSpawner(Spawner):
     def __init__(
         self,
         mode: str = "explore",
-        target: Optional[list] = None,
+        target: list | None = None,
         **_,
     ):
         self.mode = mode
diff --git a/autosampler/spawners/lof.py b/autosampler/spawners/lof.py
index 803f349..02a43a3 100644
--- a/autosampler/spawners/lof.py
+++ b/autosampler/spawners/lof.py
@@ -1,7 +1,5 @@
 from __future__ import annotations
 
-from typing import Optional
-
 import numpy as np
 
 from .base import Spawner, SpawnerFactory
@@ -20,7 +18,7 @@ def __init__(
         self,
         n_neighbors: int = 20,
         mode: str = "explore",
-        target: Optional[list] = None,
+        target: list | None = None,
         **_,
     ):
         self.n_neighbors = n_neighbors
diff --git a/autosampler/spawners/voronoi.py b/autosampler/spawners/voronoi.py
index 175ccdb..9c97f7c 100644
--- a/autosampler/spawners/voronoi.py
+++ b/autosampler/spawners/voronoi.py
@@ -5,8 +5,14 @@
 import numpy as np
 
 from autosampler.binning.spatial import VoronoiBinner
+
 from .base import SpawnerFactory
-from .density import DensitySpawner, _cumulative_points, _sample_frames, _weighted_choice
+from .density import (
+    DensitySpawner,
+    _cumulative_points,
+    _sample_frames,
+    _weighted_choice,
+)
 
 
 class VoronoiSpawner(DensitySpawner):
diff --git a/autosampler/templates.py b/autosampler/templates.py
index 33472bc..cadb4ee 100644
--- a/autosampler/templates.py
+++ b/autosampler/templates.py
@@ -66,7 +66,7 @@
 # ---- CV space: fixed physical CVs or a learned latent space ----------------
 # space_mode: fixed | pca | tica | tvae | vampnet | spib | deep-tica | deep-lda
 space_mode: vampnet
-adaptive_feature_type: distances        # distances | fitted_coords | phi_psi
+adaptive_feature_type: distances        # distances | fitted_coords | phi_psi (AIB9-only)
 retrain_freq: 5                         # retrain cadence for retrain_policy: fixed
 retrain_policy: fixed                   # fixed | vamp_adaptive (retrain on VAMP-2 drop)
 # vamp_retrain_tol: 0.1                 # relative VAMP-2 drop that triggers a retrain
@@ -135,6 +135,7 @@
 # ---- Execution: where walkers run ------------------------------------------
 execution:
   backend: local                        # local | slurm | pbs
+  # walker_timeout: 3600                 # local: kill a walker after N seconds (hang guard)
   # --- scheduler settings (slurm/pbs) ---
   # partition: gpu
   # account: my_alloc
diff --git a/autosampler/utils/math.py b/autosampler/utils/math.py
index 19de262..e8ee8f1 100644
--- a/autosampler/utils/math.py
+++ b/autosampler/utils/math.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def safe_divide(numerator: np.ndarray, denominator: np.ndarray, eps: float = 1e-8) -> np.ndarray:
     """Perform element-wise division with epsilon handling to avoid NaN/Inf."""
     denom_safe = np.where(np.abs(denominator) > eps, denominator, eps * np.sign(denominator + 1e-12))
diff --git a/autosampler/utils/seeds.py b/autosampler/utils/seeds.py
index a5aaa74..34b5736 100644
--- a/autosampler/utils/seeds.py
+++ b/autosampler/utils/seeds.py
@@ -35,6 +35,15 @@ def set_seed(self) -> None:
             torch.backends.cudnn.deterministic = True
             torch.backends.cudnn.benchmark = False
 
+        # Request deterministic algorithm implementations where available. Some
+        # CUDA kernels need this workspace setting to be deterministic; warn_only
+        # avoids hard-failing on the rare op that lacks a deterministic variant.
+        os.environ.setdefault("CUBLAS_WORKSPACE_CONFIG", ":4096:8")
+        try:
+            torch.use_deterministic_algorithms(True, warn_only=True)
+        except Exception:  # noqa: BLE001 - older torch without warn_only support
+            pass
+
         # 4. PyTorch Lightning (deep-TICA / deep-LDA), if installed.
         self._seed_lightning()
 
diff --git a/autosampler/workflows/parallel.py b/autosampler/workflows/parallel.py
index 5cd8ec5..1c789bf 100644
--- a/autosampler/workflows/parallel.py
+++ b/autosampler/workflows/parallel.py
@@ -10,7 +10,6 @@
 from __future__ import annotations
 
 from pathlib import Path
-from typing import Optional
 
 from autosampler.execution import build_walker_tasks, make_backend
 
@@ -25,7 +24,7 @@ def run_iteration_parallel(
     outdir: Path,
     iteration: int,
     max_workers: int = 8,
-    gpu_ids: Optional[list[int]] = None,
+    gpu_ids: list[int] | None = None,
     execution=None,
 ) -> list:
     """Execute walker production runs via the configured execution backend."""
diff --git a/docs/api.md b/docs/api.md
new file mode 100644
index 0000000..9bed7b3
--- /dev/null
+++ b/docs/api.md
@@ -0,0 +1,28 @@
+# API reference
+
+Auto-generated from the source docstrings. These are the modules a user or
+referee is most likely to call directly; the full package has more.
+
+## Configuration
+
+::: autosampler.config.AutoSamplerConfig
+
+## MSM estimation & convergence
+
+::: autosampler.msm.estimator
+
+::: autosampler.msm.convergence
+
+## Collective-variable feature selection
+
+::: autosampler.spaces.feature_selection
+
+## Adaptive binning
+
+::: autosampler.binning.adaptive
+
+## Analysis & plotting
+
+::: autosampler.analysis.data
+
+::: autosampler.analysis.plots
diff --git a/docs/cli.md b/docs/cli.md
new file mode 100644
index 0000000..a8cd5f2
--- /dev/null
+++ b/docs/cli.md
@@ -0,0 +1,79 @@
+# CLI reference
+
+AutoSampler installs six console commands (see `pyproject.toml [project.scripts]`).
+
+## `autosampler` (alias `autosampler-run`)
+
+Run the adaptive sampling loop from a config file.
+
+```bash
+autosampler --config CONFIG.yaml [--iterations N] [--resume latest|N] [--check] [--log-level LEVEL]
+```
+
+| Flag | Default | Description |
+| --- | --- | --- |
+| `--config` | *(required)* | Path to the YAML input file. |
+| `--iterations` | `1` | Number of iterations to run (`>= 0`). |
+| `--resume` | — | Resume from the latest checkpoint (`latest`) or a specific iteration (`N`). |
+| `--check` | off | Preflight only: validate inputs/engine/settings, then exit without running MD. |
+| `--log-level` | `INFO` | Python logging level (e.g. `DEBUG`, `WARNING`). |
+
+## `autosampler-init`
+
+Write a fully-annotated starter input file.
+
+```bash
+autosampler-init [-o OUTPUT] [--force]
+```
+
+| Flag | Default | Description |
+| --- | --- | --- |
+| `-o`, `--output` | `config.yaml` | Where to write the template. |
+| `--force` | off | Overwrite an existing file. |
+
+## `autosampler-analyze`
+
+Produce a multi-panel MSM convergence report from a run directory.
+
+```bash
+autosampler-analyze --run-dir RUN_DIR [--outfile FILE] [--temperature K]
+```
+
+| Flag | Default | Description |
+| --- | --- | --- |
+| `--run-dir` | *(required)* | Run output directory (contains `iter_*/msm.npz`). |
+| `--outfile` | — | Output image path for the report. |
+| `--temperature` | `300.0` | Temperature (K) for free-energy conversion. |
+
+## `autosampler-log`
+
+Write/extend an exploration log (per-iteration CV-bin occupancy) for a run.
+
+```bash
+autosampler-log --run-dir RUN_DIR [--config CONFIG] [--output FILE] \
+    [--n-bins ...] [--min-values ...] [--max-values ...] [--append]
+```
+
+| Flag | Default | Description |
+| --- | --- | --- |
+| `--run-dir` | *(required)* | Run output directory. |
+| `--config` | — | Config to read binning bounds from. |
+| `--output` | — | Output log path. |
+| `--n-bins` / `--min-values` / `--max-values` | — | Override the CV grid (comma-separated lists). |
+| `--append` | off | Append to an existing log. |
+
+## `autosampler-path`
+
+Reconstruct a connected trajectory between two CV points from a run's frame
+lineage.
+
+```bash
+autosampler-path --run-dir RUN_DIR --topology TOP --start "x,y" --end "x,y" [--output OUT]
+```
+
+| Flag | Default | Description |
+| --- | --- | --- |
+| `--run-dir` | *(required)* | Run output directory. |
+| `--topology` | *(required)* | Topology for writing the connected trajectory. |
+| `--start` / `--end` | — | CV coordinates (comma-separated) of the path endpoints. |
+| `--output` | — | Output trajectory path. |
diff --git a/docs/configuration.md b/docs/configuration.md
index 66e0617..9664bd4 100644
--- a/docs/configuration.md
+++ b/docs/configuration.md
@@ -44,8 +44,21 @@ at startup. Below, only non-obvious defaults are noted — see
 | `walker` / `step` / `stride` | `10` / `10000` / `100` | Walkers per iteration / MD steps / save interval. |
 | `max_workers` | `4` | Concurrent walkers (local backend). |
 | `voronoi_clusters` | `150` | Cells / microstates (also used by the MSM spawner). |
+| `target` | — | CV-space target `[x, y, …]` when `search_mode: target`. |
+| `recent_density_window` | `5` | Bins sampled in the last N iterations are down-weighted (`density`). |
+| `lof_neighbors` | `20` | Neighbours for the LOF spawner. |
+| `voronoi_periodic` | `false` | Wrap Voronoi cells periodically. |
+| `voronoi_grid_size` | `250` | Grid resolution for Voronoi cell-area estimation. |
+| `voronoi_max_clusters` | `5000` | Upper bound on auto-grown Voronoi cells. |
+| `resolution_check_patience` | `5` | Iterations of bin-occupancy stall before refining the grid. |
+| `resolution_max_bins` | `150` | Upper bound on per-axis bins when auto-refining. |
 | `convergence_patience` | `0` | Bin-occupancy stall patience (legacy convergence). |
 
+!!! note "Walker timeout (local backend)"
+    `execution.walker_timeout` (seconds) kills a walker that runs longer than the
+    limit and marks the batch failed — a guard against a hung in-process OpenMM
+    walker. Off by default.
+
 ## `space_mode` and adaptive model
 
 `space_mode`: `fixed` \| `pca` \| `tica` \| `tvae` \| `vampnet` \| `spib` \|
@@ -53,7 +66,7 @@ at startup. Below, only non-obvious defaults are noted — see
 
 | Key | Default | Description |
 | --- | --- | --- |
-| `adaptive_feature_type` | `distances` | `distances` \| `fitted_coords` \| `phi_psi`. |
+| `adaptive_feature_type` | `distances` | `distances` \| `fitted_coords` \| `phi_psi`. **Note:** `phi_psi` is currently specific to the AIB9 peptide (it expects 9 `resname AIB` residues) and will raise on other systems — use `distances` or `fitted_coords` for general systems. |
 | `retrain_freq` | `1` | Retrain the CV every N iterations (`fixed` policy). |
 | `retrain_policy` | `fixed` | `fixed` or `vamp_adaptive` (retrain on VAMP-2 drop). |
 | `vamp_retrain_tol` | `0.1` | Relative VAMP-2 drop that triggers a retrain. |
diff --git a/docs/cv_methods.md b/docs/cv_methods.md
index 6f0be94..0813037 100644
--- a/docs/cv_methods.md
+++ b/docs/cv_methods.md
@@ -1,5 +1,8 @@
 # Collective variables
 
+> Primary references for TICA, VAMPNets, SPIB, deep-TICA, etc. are collected on
+> the [References](references.md) page.
+
 AutoSampler can sample in **fixed** physical CVs or **learn** CVs on the fly.
 The available learned methods live in a single registry
 (`autosampler/spaces/registry.py`), which also tracks each method's backend and
diff --git a/docs/notebook_tutorial.md b/docs/notebook_tutorial.md
index c16bcd8..1ceabe8 100644
--- a/docs/notebook_tutorial.md
+++ b/docs/notebook_tutorial.md
@@ -1,7 +1,7 @@
 # Notebook tutorial
 
 A runnable Jupyter notebook with **rendered plots** lives at
-[`examples/notebooks/adaptive_msm_tutorial.ipynb`](https://github.com/TeamSuman/AutoSampler/blob/devel/examples/notebooks/adaptive_msm_tutorial.ipynb).
+[`examples/notebooks/adaptive_msm_tutorial.ipynb`](https://github.com/TeamSuman/AutoSampler/blob/main/examples/notebooks/adaptive_msm_tutorial.ipynb).
 It uses small synthetic examples so every figure renders in seconds without
 running molecular dynamics.
 
diff --git a/docs/references.md b/docs/references.md
new file mode 100644
index 0000000..d22687f
--- /dev/null
+++ b/docs/references.md
@@ -0,0 +1,50 @@
+# References
+
+The methods AutoSampler builds on, with primary references. Cite the relevant
+ones alongside AutoSampler (see [`CITATION.cff`](https://github.com/TeamSuman/AutoSampler/blob/main/CITATION.cff)).
+
+## Markov state models
+
+- Prinz, J.-H. et al. *Markov models of molecular kinetics: Generation and
+  validation.* J. Chem. Phys. **134**, 174105 (2011).
+- Husic, B. E. & Pande, V. S. *Markov state models: From an art to a science.*
+  J. Am. Chem. Soc. **140**, 2386–2396 (2018).
+- Röblitz, S. & Weber, M. *Fuzzy spectral clustering by PCCA+.* Adv. Data Anal.
+  Classif. **7**, 147–179 (2013). (PCCA+ metastable decomposition.)
+
+## Dimensionality reduction / collective variables
+
+- Pérez-Hernández, G. et al. *Identification of slow molecular order parameters
+  for Markov model construction (TICA).* J. Chem. Phys. **139**, 015102 (2013).
+- Schwantes, C. R. & Pande, V. S. *Improvements in Markov state model
+  construction reveal many non-native interactions in the folding of NTL9.*
+  J. Chem. Theory Comput. **9**, 2000–2009 (2013). (TICA.)
+- Wu, H. & Noé, F. *Variational approach for learning Markov processes from time
+  series data (VAMP / VAMP-2 score).* J. Nonlinear Sci. **30**, 23–66 (2020).
+- Mardt, A. et al. *VAMPnets for deep learning of molecular kinetics.* Nat.
+  Commun. **9**, 5 (2018). (VAMPNet.)
+- Bonati, L., Piccini, G. & Parrinello, M. *Deep learning the slow modes for rare
+  events sampling (Deep-TICA).* PNAS **118**, e2113533118 (2021).
+- Wang, D. & Tiwary, P. *State predictive information bottleneck (SPIB).*
+  J. Chem. Phys. **154**, 134111 (2021).
+- Time-lagged (variational) autoencoders: Wehmeyer, C. & Noé, F. *Time-lagged
+  autoencoders.* J. Chem. Phys. **148**, 241703 (2018).
+
+## Adaptive sampling & weighted ensemble
+
+- Huber, G. A. & Kim, S. *Weighted-ensemble Brownian dynamics simulations.*
+  Biophys. J. **70**, 97–110 (1996).
+- Zuckerman, D. M. & Chong, L. T. *Weighted ensemble simulation: Review of
+  methodology, applications, and software.* Annu. Rev. Biophys. **46**, 43–57
+  (2017).
+- Minimal adaptive binning (MAB) for weighted ensemble: Torrillo, P. A.,
+  DeGrave, A. J. & Bhatt, J. *A minimal, adaptive binning scheme for weighted
+  ensemble simulations.* J. Phys. Chem. A **125**, 1642–1649 (2021).
+
+## Software
+
+- Hoffmann, M. et al. *Deeptime: a Python library for machine learning dynamical
+  models from time series data.* Mach. Learn.: Sci. Technol. **3**, 015009 (2022).
+- Eastman, P. et al. *OpenMM 8.* J. Phys. Chem. B **128**, 109–116 (2024).
+- Michaud-Agrawal, N. et al. *MDAnalysis: A toolkit for the analysis of molecular
+  dynamics simulations.* J. Comput. Chem. **32**, 2319–2327 (2011).
diff --git a/env.yml b/env.yml
index cc7b811..1e5852e 100644
--- a/env.yml
+++ b/env.yml
@@ -8,7 +8,8 @@ dependencies:
   - wheel
   - numpy>=1.23
   - scipy>=1.9
-  - pydantic>=1.10
+  - shapely>=2.0          # required by the Voronoi spawner/binner (core dependency)
+  - pydantic>=2.0         # code uses Pydantic v2 (field_validator/model_validator)
   - pyyaml>=6.0
   - scikit-learn>=1.2
   - mdanalysis>=2.5
@@ -18,4 +19,7 @@ dependencies:
   - matplotlib>=3.7
   - pytest>=7.0
   - pip:
+      # Optional deep-CV backends (space_mode: deep-tica / deep-lda).
+      - lightning>=2.0
+      - mlcolvar>=1.0
       - -e .
diff --git a/examples/AIB9/config_deep_tica.yaml b/examples/AIB9/config_deep_tica.yaml
new file mode 100644
index 0000000..5611b00
--- /dev/null
+++ b/examples/AIB9/config_deep_tica.yaml
@@ -0,0 +1,41 @@
+# deep-TICA collective variable on AIB9 (requires: pip install '.[deep-tica]').
+system:
+  conf_file: aib9_equilibrated.pdb
+  top_file: aib9_equilibrated.pdb
+  topology: charmm
+  system_file: system.py
+  trajectory_topology_file: aib9_equilibrated.pdb
+  feature_selection: "resname AIB and name N CA C O CB1 CB2"
+
+engine:
+  md_engine: openmm
+  platform_name: CUDA          # use CPU if no GPU (slower)
+  precision: mixed
+  npt: true
+  temperature: 400.0
+  dt: 0.002
+
+spawning:
+  spawn_scheme: density
+  search_mode: explore
+  walker: 10
+  step: 5000
+  stride: 50
+  max_workers: 1
+
+space_mode: deep-tica
+adaptive_feature_type: distances
+retrain_freq: 5
+aggregate_memory: true
+max_adaptive_memory_frames: 5000
+adaptive_model:
+  lagtime: 5
+  latent_dim: 2
+  epochs: 50
+  deep_tica_hidden_dims: [64, 32]
+
+outdir: runs/aib9_deep_tica
+random_seed: 7
+checkpoint_freq: 1
+save_features: true
+n_bins: [30, 30]
diff --git a/examples/AIB9/config_pbs.yaml b/examples/AIB9/config_pbs.yaml
new file mode 100644
index 0000000..9254485
--- /dev/null
+++ b/examples/AIB9/config_pbs.yaml
@@ -0,0 +1,43 @@
+# PBS/Torque HPC execution: one array job per iteration. See docs/execution.md.
+system:
+  conf_file: aib9_equilibrated.pdb
+  top_file: aib9_equilibrated.pdb
+  topology: charmm
+  system_file: system.py
+  trajectory_topology_file: aib9_equilibrated.pdb
+  feature_selection: "resname AIB and name N CA C O CB1 CB2"
+
+engine:
+  md_engine: openmm
+  platform_name: CUDA          # use CPU if no GPU (slower)
+  precision: mixed
+  npt: true
+  temperature: 400.0
+  dt: 0.002
+
+spawning:
+  spawn_scheme: density
+  search_mode: explore
+  walker: 32
+  step: 5000
+  stride: 50
+
+space_mode: tica
+adaptive_feature_type: distances
+retrain_freq: 5
+
+execution:
+  backend: pbs
+  partition: gpu               # PBS queue
+  walltime: "02:00:00"
+  cpus_per_task: 8
+  gpus_per_task: 1
+  memory: "16G"
+  max_retries: 2
+  module_loads:
+    - "module load cuda/12.2"
+
+outdir: runs/aib9_pbs
+random_seed: 7
+checkpoint_freq: 1
+n_bins: [30, 30]
diff --git a/examples/AIB9/config_spib.yaml b/examples/AIB9/config_spib.yaml
new file mode 100644
index 0000000..72af9d5
--- /dev/null
+++ b/examples/AIB9/config_spib.yaml
@@ -0,0 +1,42 @@
+# SPIB (State Predictive Information Bottleneck) collective variable on AIB9.
+system:
+  conf_file: aib9_equilibrated.pdb
+  top_file: aib9_equilibrated.pdb
+  topology: charmm
+  system_file: system.py
+  trajectory_topology_file: aib9_equilibrated.pdb
+  feature_selection: "resname AIB and name N CA C O CB1 CB2"
+
+engine:
+  md_engine: openmm
+  platform_name: CUDA          # use CPU if no GPU (slower)
+  precision: mixed
+  npt: true
+  temperature: 400.0
+  dt: 0.002
+
+spawning:
+  spawn_scheme: density
+  search_mode: explore
+  walker: 10
+  step: 5000
+  stride: 50
+  max_workers: 1
+
+space_mode: spib
+adaptive_feature_type: distances
+retrain_freq: 5
+aggregate_memory: true
+max_adaptive_memory_frames: 5000
+adaptive_model:
+  lagtime: 5
+  latent_dim: 2
+  epochs: 50
+  spib_n_states: 10
+  spib_beta: 0.001
+
+outdir: runs/aib9_spib
+random_seed: 7
+checkpoint_freq: 1
+save_features: true
+n_bins: [30, 30]
diff --git a/examples/AIB9/config_target.yaml b/examples/AIB9/config_target.yaml
new file mode 100644
index 0000000..341d035
--- /dev/null
+++ b/examples/AIB9/config_target.yaml
@@ -0,0 +1,37 @@
+# Target-mode search: drive sampling toward a chosen point in the learned CV space.
+system:
+  conf_file: aib9_equilibrated.pdb
+  top_file: aib9_equilibrated.pdb
+  topology: charmm
+  system_file: system.py
+  trajectory_topology_file: aib9_equilibrated.pdb
+  feature_selection: "resname AIB and name N CA C O CB1 CB2"
+
+engine:
+  md_engine: openmm
+  platform_name: CUDA          # use CPU if no GPU (slower)
+  precision: mixed
+  npt: true
+  temperature: 400.0
+  dt: 0.002
+
+spawning:
+  spawn_scheme: density
+  spawn_type: probabilistic
+  search_mode: target          # bias spawning toward `target`
+  target: [0.0, 0.0]           # target coordinates in the 2-D latent CV space
+  walker: 10
+  step: 5000
+  stride: 50
+  max_workers: 1
+
+space_mode: tica
+adaptive_feature_type: distances
+retrain_freq: 5
+aggregate_memory: true
+max_adaptive_memory_frames: 5000
+
+outdir: runs/aib9_target
+random_seed: 7
+checkpoint_freq: 1
+n_bins: [30, 30]
diff --git a/examples/AIB9/config_we.yaml b/examples/AIB9/config_we.yaml
new file mode 100644
index 0000000..0fcc3b2
--- /dev/null
+++ b/examples/AIB9/config_we.yaml
@@ -0,0 +1,36 @@
+# Weighted-ensemble (WE) resampling spawner on AIB9 (fixed phi/psi CVs).
+system:
+  conf_file: aib9_equilibrated.pdb
+  top_file: aib9_equilibrated.pdb
+  topology: charmm
+  system_file: system.py
+  trajectory_topology_file: aib9_equilibrated.pdb
+  feature_selection: "resname AIB and name N CA C O CB1 CB2"
+
+engine:
+  md_engine: openmm
+  platform_name: CUDA          # use CPU if no GPU (slower)
+  precision: mixed
+  npt: true
+  temperature: 400.0
+  dt: 0.002
+
+spawning:
+  spawn_scheme: we             # weight-conserving split/merge per occupied bin
+  we_target_per_bin: 4
+  search_mode: explore
+  walker: 16
+  step: 5000
+  stride: 50
+  max_workers: 1
+
+space_mode: tica
+adaptive_feature_type: distances
+retrain_freq: 5
+aggregate_memory: true
+max_adaptive_memory_frames: 5000
+
+outdir: runs/aib9_we
+random_seed: 7
+checkpoint_freq: 1
+n_bins: [30, 30]
diff --git a/examples/AlaD/config.yaml b/examples/AlaD/config.yaml
index 6c7caa2..5778225 100644
--- a/examples/AlaD/config.yaml
+++ b/examples/AlaD/config.yaml
@@ -9,9 +9,13 @@ system:
 engine:
   # md_engine: gromacs
   md_engine: openmm
-  platform: CUDA
+  platform_name: CUDA            # CUDA | CPU | OpenCL | Reference  (use CPU if no GPU)
   gromacs_executable: gmx
-  gromacs_include_dir: /home/dm/Soft/GMX26/share/gromacs/top
+  # OpenMM parses this GROMACS topology, so it needs the GROMACS force-field
+  # directory (amber99sb.ff/...) to resolve the #include lines. Point this at
+  # your install's `share/gromacs/top` (or run `gmx pdb2gmx`/`grompp -pp` to make
+  # a self-contained processed topology).
+  gromacs_include_dir: /path/to/gromacs/share/gromacs/top
   temperature: 300.0
   pressure: 1.0
   dt: 0.002
diff --git a/examples/AlaD/config_voronoi.yaml b/examples/AlaD/config_voronoi.yaml
index c00817d..d234e3e 100644
--- a/examples/AlaD/config_voronoi.yaml
+++ b/examples/AlaD/config_voronoi.yaml
@@ -9,9 +9,11 @@ system:
 engine:
   # md_engine: gromacs
   md_engine: openmm
-  platform: CUDA
+  platform_name: CUDA            # CUDA | CPU | OpenCL | Reference  (use CPU if no GPU)
   gromacs_executable: gmx
-  gromacs_include_dir: /home/dm/Soft/GMX26/share/gromacs/top
+  # Point this at your GROMACS install's `share/gromacs/top` so OpenMM can resolve
+  # the amber99sb.ff #include lines in topol.top.
+  gromacs_include_dir: /path/to/gromacs/share/gromacs/top
   temperature: 300.0
   pressure: 1.0
   dt: 0.002
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..04f8478
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,62 @@
+# AutoSampler examples
+
+Each example is a self-contained config plus the assets it needs. Paths inside a
+config are resolved relative to that config file. Validate any example without
+running MD using `--check`:
+
+```bash
+autosampler --config examples/<path>/config.yaml --check
+```
+
+## Start here (laptop-friendly, no GPU)
+
+| Example | Demonstrates | Requirements |
+| --- | --- | --- |
+| [`alanine_dipeptide/config.yaml`](alanine_dipeptide/) | **Hello world** — fixed phi/psi CVs, density spawning, the full adaptive loop end-to-end | CPU only; vacuum Amber14 system, **no external force-field files**. Runs in minutes. |
+
+The alanine-dipeptide system (`system.xml`) is built from the bundled 22-atom
+`structure.pdb` (the canonical OpenMM test structure, MIT) by
+`build_system.py` — rerun it to regenerate the asset. This is the recommended
+first run:
+
+```bash
+autosampler --config examples/alanine_dipeptide/config.yaml --iterations 5
+```
+
+## AIB9 peptide (uses the bundled `aib9_system.xml`)
+
+| Example | Demonstrates | Requirements |
+| --- | --- | --- |
+| `AIB9/config_fixed_phi_psi.yaml` | Fixed CVs on AIB9 | GPU recommended (solvated system) |
+| `AIB9/config_phi_psi.yaml` | phi/psi projection | GPU recommended |
+| `AIB9/config_adaptive.yaml` | Learned CV (TVAE) | GPU recommended |
+| `AIB9/config_msm_vampnet.yaml` | VAMPNet CV + MSM convergence + MSM spawning + gradient binning | GPU recommended; `pip install '.[deep-tica]'` |
+| `AIB9/config_msm_feature_selection.yaml` | VAMP-2 feature selection + MSM | GPU recommended |
+| `AIB9/config_spib.yaml` | SPIB learned CV | GPU recommended |
+| `AIB9/config_deep_tica.yaml` | deep-TICA learned CV | GPU recommended; `pip install '.[deep-tica]'` |
+| `AIB9/config_we.yaml` | Weighted-ensemble (WE) spawner | GPU recommended |
+| `AIB9/config_target.yaml` | Target-mode search toward a CV point | GPU recommended |
+| `AIB9/config_pbs.yaml` | PBS/Torque HPC backend (array job per iteration) | A PBS cluster |
+| `AIB9/config_tda_phi_sweep.yaml` | TDA phi sweep | Needs external Zenodo assets (see the config header) |
+
+## Alanine dipeptide (GROMACS topology via OpenMM)
+
+| Example | Demonstrates | Requirements |
+| --- | --- | --- |
+| `AlaD/config.yaml` | Fixed phi/psi, density spawner | A GROMACS install (its `share/gromacs/top` for the amber99sb.ff includes); set `engine.gromacs_include_dir`. GPU recommended. |
+| `AlaD/config_voronoi.yaml` | Voronoi spawner | Same as above |
+
+## Run scripts
+
+- `run_local.sh` — run on a workstation (local backend).
+- `slurm_submit.sh` / `pbs_submit.sh` — submit on an HPC cluster (see
+  [docs/execution.md](../docs/execution.md)).
+
+## Notebook
+
+- `notebooks/adaptive_msm_tutorial.ipynb` — a rendered, GPU-free walkthrough of
+  the components (feature selection → MSM → convergence → weighted ensemble →
+  analysis) on synthetic data.
+
+> Run outputs are written under each example's `outdir` (e.g. `runs/...`) and are
+> not committed.
diff --git a/examples/alanine_dipeptide/build_system.py b/examples/alanine_dipeptide/build_system.py
new file mode 100644
index 0000000..5455776
--- /dev/null
+++ b/examples/alanine_dipeptide/build_system.py
@@ -0,0 +1,35 @@
+"""Generate the vacuum alanine-dipeptide OpenMM system for the hello-world example.
+
+Builds an Amber14 vacuum ``System`` from ``structure.pdb`` (Ac-Ala-NMe, 22 atoms)
+and serialises it to ``system.xml``. Run once to (re)generate the committed asset:
+
+    python examples/alanine_dipeptide/build_system.py
+
+``structure.pdb`` is the canonical alanine-dipeptide test structure from the
+OpenMM project (MIT licensed). Only OpenMM is required — no external force-field
+files, no GPU.
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+from openmm import XmlSerializer
+from openmm.app import ForceField, NoCutoff, PDBFile
+
+BASE_DIR = Path(__file__).resolve().parent
+
+
+def main() -> None:
+    pdb = PDBFile(str(BASE_DIR / "structure.pdb"))
+    forcefield = ForceField("amber14-all.xml")  # vacuum: no water model
+    system = forcefield.createSystem(
+        pdb.topology, nonbondedMethod=NoCutoff, constraints=None
+    )
+    out = BASE_DIR / "system.xml"
+    out.write_text(XmlSerializer.serialize(system))
+    print("WROTE", out, "with", system.getNumParticles(), "particles")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/alanine_dipeptide/config.yaml b/examples/alanine_dipeptide/config.yaml
new file mode 100644
index 0000000..2064f99
--- /dev/null
+++ b/examples/alanine_dipeptide/config.yaml
@@ -0,0 +1,44 @@
+# Alanine-dipeptide hello-world: a self-contained, CPU-only, GPU-free example.
+#
+# Uses a vacuum Amber14 system built from a 22-atom structure (no external
+# force-field files, no solvent). Runs in minutes on a laptop. CV space is the
+# phi/psi dihedral pair; walkers are spawned by density over a phi/psi grid.
+#
+# Regenerate system.xml with:  python examples/alanine_dipeptide/build_system.py
+# Validate:                    autosampler --config examples/alanine_dipeptide/config.yaml --check
+# Run:                         autosampler --config examples/alanine_dipeptide/config.yaml --iterations 5
+
+system:
+  conf_file: structure.pdb
+  top_file: structure.pdb
+  topology: amber
+  system_file: system.py            # builds the OpenMM System + integrator
+  project_file: project_phi_psi.py  # defines extract_cvs (phi, psi)
+  trajectory_topology_file: structure.pdb
+  feature_selection: "protein and not (type H)"
+
+engine:
+  md_engine: openmm
+  platform_name: CPU                # laptop-friendly; no GPU required
+  temperature: 300.0
+  dt: 0.002
+  npt: false
+
+spawning:
+  spawn_scheme: density
+  spawn_type: hard
+  search_mode: explore
+  walker: 4                         # small batch
+  step: 1000                        # 2 ps per walker
+  stride: 100
+  max_workers: 2
+
+space_mode: fixed
+n_bins: [18, 18]
+min_values: [-3.141592653589793, -3.141592653589793]
+max_values: [3.141592653589793, 3.141592653589793]
+
+outdir: runs/alanine_dipeptide_hello
+random_seed: 42
+checkpoint_freq: 1
+save_features: false
diff --git a/examples/alanine_dipeptide/project_phi_psi.py b/examples/alanine_dipeptide/project_phi_psi.py
new file mode 100644
index 0000000..330c6b2
--- /dev/null
+++ b/examples/alanine_dipeptide/project_phi_psi.py
@@ -0,0 +1,46 @@
+"""Alanine dipeptide phi/psi CV extractor for AutoSampler."""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def extract_cvs(trajectories, top_file, conf_file):
+    import MDAnalysis as mda
+    from MDAnalysis.lib.distances import calc_dihedrals
+
+    u = mda.Universe(conf_file, trajectories)
+    try:
+        ace_c = _single_atom(u, "resname ACE and name C")
+        ala_n = _single_atom(u, "resname ALA and name N")
+        ala_ca = _single_atom(u, "resname ALA and name CA")
+        ala_c = _single_atom(u, "resname ALA and name C")
+        nme_n = _single_atom(u, "resname NME and name N")
+
+        cvs = np.zeros((u.trajectory.n_frames, 2), dtype=np.float32)
+        for frame_index, ts in enumerate(u.trajectory):
+            phi = calc_dihedrals(
+                ace_c.positions,
+                ala_n.positions,
+                ala_ca.positions,
+                ala_c.positions,
+                box=ts.dimensions,
+            )[0]
+            psi = calc_dihedrals(
+                ala_n.positions,
+                ala_ca.positions,
+                ala_c.positions,
+                nme_n.positions,
+                box=ts.dimensions,
+            )[0]
+            cvs[frame_index] = [phi, psi]
+        return cvs
+    finally:
+        u.trajectory.close()
+
+
+def _single_atom(universe, selection):
+    atoms = universe.select_atoms(selection)
+    if atoms.n_atoms != 1:
+        raise ValueError(f"Selection {selection!r} matched {atoms.n_atoms} atoms.")
+    return atoms
diff --git a/examples/alanine_dipeptide/structure.pdb b/examples/alanine_dipeptide/structure.pdb
new file mode 100644
index 0000000..32ba37b
--- /dev/null
+++ b/examples/alanine_dipeptide/structure.pdb
@@ -0,0 +1,25 @@
+REMARK  ACE                                                         
+ATOM      1 1HH3 ACE     1       2.000   1.000  -0.000
+ATOM      2  CH3 ACE     1       2.000   2.090   0.000
+ATOM      3 2HH3 ACE     1       1.486   2.454   0.890
+ATOM      4 3HH3 ACE     1       1.486   2.454  -0.890
+ATOM      5  C   ACE     1       3.427   2.641  -0.000
+ATOM      6  O   ACE     1       4.391   1.877  -0.000
+ATOM      7  N   ALA     2       3.555   3.970  -0.000
+ATOM      8  H   ALA     2       2.733   4.556  -0.000
+ATOM      9  CA  ALA     2       4.853   4.614  -0.000
+ATOM     10  HA  ALA     2       5.408   4.316   0.890
+ATOM     11  CB  ALA     2       5.661   4.221  -1.232
+ATOM     12 1HB  ALA     2       5.123   4.521  -2.131
+ATOM     13 2HB  ALA     2       6.630   4.719  -1.206
+ATOM     14 3HB  ALA     2       5.809   3.141  -1.241
+ATOM     15  C   ALA     2       4.713   6.129   0.000
+ATOM     16  O   ALA     2       3.601   6.653   0.000
+ATOM     17  N   NME     3       5.846   6.835   0.000
+ATOM     18  H   NME     3       6.737   6.359  -0.000
+ATOM     19  CH3 NME     3       5.846   8.284   0.000
+ATOM     20 1HH3 NME     3       4.819   8.648   0.000
+ATOM     21 2HH3 NME     3       6.360   8.648   0.890
+ATOM     22 3HH3 NME     3       6.360   8.648  -0.890
+TER   
+END   
diff --git a/examples/alanine_dipeptide/system.py b/examples/alanine_dipeptide/system.py
new file mode 100644
index 0000000..6b403da
--- /dev/null
+++ b/examples/alanine_dipeptide/system.py
@@ -0,0 +1,22 @@
+"""OpenMM system loader for the alanine-dipeptide hello-world example."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+from openmm import LangevinMiddleIntegrator, XmlSerializer
+from openmm.unit import kelvin, picosecond, picoseconds
+
+BASE_DIR = Path(__file__).resolve().parent
+
+
+def make_system(_topology_source, temp=300.0, dt=0.002):
+    """Load the vacuum alanine-dipeptide system built by build_system.py."""
+    with (BASE_DIR / "system.xml").open() as handle:
+        system = XmlSerializer.deserialize(handle.read())
+    integrator = LangevinMiddleIntegrator(
+        temp * kelvin,
+        1.0 / picosecond,
+        dt * picoseconds,
+    )
+    return system, integrator
diff --git a/examples/alanine_dipeptide/system.xml b/examples/alanine_dipeptide/system.xml
new file mode 100644
index 0000000..5386e36
--- /dev/null
+++ b/examples/alanine_dipeptide/system.xml
@@ -0,0 +1,276 @@
+<?xml version="1.0" ?>
+<System openmmVersion="8.5.2" type="System" version="1">
+	<PeriodicBoxVectors>
+		<A x="2" y="0" z="0"/>
+		<B x="0" y="2" z="0"/>
+		<C x="0" y="0" z="2"/>
+	</PeriodicBoxVectors>
+	<Particles>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="16"/>
+		<Particle mass="14.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="1.008"/>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="16"/>
+		<Particle mass="14.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="12.01"/>
+		<Particle mass="1.008"/>
+		<Particle mass="1.008"/>
+		<Particle mass="1.008"/>
+	</Particles>
+	<Constraints/>
+	<Forces>
+		<Force forceGroup="0" name="HarmonicBondForce" type="HarmonicBondForce" usesPeriodic="0" version="2">
+			<Bonds>
+				<Bond d=".1522" k="265265.6" p1="4" p2="1"/>
+				<Bond d=".12290000000000001" k="476975.9999999999" p1="4" p2="5"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="1" p2="0"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="1" p2="2"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="1" p2="3"/>
+				<Bond d=".1335" k="410031.99999999994" p1="4" p2="6"/>
+				<Bond d=".1522" k="265265.6" p1="14" p2="8"/>
+				<Bond d=".12290000000000001" k="476975.9999999999" p1="14" p2="15"/>
+				<Bond d=".1526" k="259407.99999999994" p1="8" p2="10"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="8" p2="9"/>
+				<Bond d=".1449" k="282001.5999999999" p1="8" p2="6"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="10" p2="11"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="10" p2="12"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="10" p2="13"/>
+				<Bond d=".101" k="363171.19999999995" p1="7" p2="6"/>
+				<Bond d=".1335" k="410031.99999999994" p1="14" p2="16"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="18" p2="19"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="18" p2="20"/>
+				<Bond d=".10900000000000001" k="284511.99999999994" p1="18" p2="21"/>
+				<Bond d=".1449" k="282001.5999999999" p1="18" p2="16"/>
+				<Bond d=".101" k="363171.19999999995" p1="17" p2="16"/>
+			</Bonds>
+		</Force>
+		<Force forceGroup="0" name="PeriodicTorsionForce" type="PeriodicTorsionForce" usesPeriodic="0" version="2">
+			<Torsions>
+				<Torsion k="3.3472000000000004" p1="0" p2="1" p3="4" p4="5" periodicity="1" phase="0"/>
+				<Torsion k=".33472" p1="0" p2="1" p3="4" p4="5" periodicity="3" phase="3.141592653589793"/>
+				<Torsion k="10.46" p1="1" p2="4" p3="6" p4="7" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="10.46" p1="1" p2="4" p3="6" p4="8" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="3.3472000000000004" p1="2" p2="1" p3="4" p4="5" periodicity="1" phase="0"/>
+				<Torsion k=".33472" p1="2" p2="1" p3="4" p4="5" periodicity="3" phase="3.141592653589793"/>
+				<Torsion k="3.3472000000000004" p1="3" p2="1" p3="4" p4="5" periodicity="1" phase="0"/>
+				<Torsion k=".33472" p1="3" p2="1" p3="4" p4="5" periodicity="3" phase="3.141592653589793"/>
+				<Torsion k="3.3472000000000004" p1="4" p2="6" p3="8" p4="10" periodicity="3" phase="0"/>
+				<Torsion k="7.5312" p1="4" p2="6" p3="8" p4="10" periodicity="2" phase="0"/>
+				<Torsion k="8.368" p1="4" p2="6" p3="8" p4="10" periodicity="1" phase="0"/>
+				<Torsion k="1.75728" p1="4" p2="6" p3="8" p4="14" periodicity="3" phase="0"/>
+				<Torsion k="1.12968" p1="4" p2="6" p3="8" p4="14" periodicity="2" phase="0"/>
+				<Torsion k="10.46" p1="5" p2="4" p3="6" p4="7" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="8.368" p1="5" p2="4" p3="6" p4="7" periodicity="1" phase="0"/>
+				<Torsion k="10.46" p1="5" p2="4" p3="6" p4="8" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k=".6508444444444444" p1="6" p2="8" p3="10" p4="11" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="6" p2="8" p3="10" p4="12" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="6" p2="8" p3="10" p4="13" periodicity="3" phase="0"/>
+				<Torsion k="2.3012" p1="6" p2="8" p3="14" p4="16" periodicity="3" phase="3.141592653589793"/>
+				<Torsion k="6.610720000000001" p1="6" p2="8" p3="14" p4="16" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="1.8828" p1="6" p2="8" p3="14" p4="16" periodicity="1" phase="3.141592653589793"/>
+				<Torsion k="10.46" p1="8" p2="14" p3="16" p4="17" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="10.46" p1="8" p2="14" p3="16" p4="18" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k=".6508444444444444" p1="9" p2="8" p3="10" p4="11" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="9" p2="8" p3="10" p4="12" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="9" p2="8" p3="10" p4="13" periodicity="3" phase="0"/>
+				<Torsion k="3.3472000000000004" p1="9" p2="8" p3="14" p4="15" periodicity="1" phase="0"/>
+				<Torsion k=".33472" p1="9" p2="8" p3="14" p4="15" periodicity="3" phase="3.141592653589793"/>
+				<Torsion k="1.6736000000000002" p1="10" p2="8" p3="14" p4="16" periodicity="3" phase="0"/>
+				<Torsion k=".8368000000000001" p1="10" p2="8" p3="14" p4="16" periodicity="2" phase="0"/>
+				<Torsion k=".8368000000000001" p1="10" p2="8" p3="14" p4="16" periodicity="1" phase="0"/>
+				<Torsion k=".6508444444444444" p1="11" p2="10" p3="8" p4="14" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="12" p2="10" p3="8" p4="14" periodicity="3" phase="0"/>
+				<Torsion k=".6508444444444444" p1="13" p2="10" p3="8" p4="14" periodicity="3" phase="0"/>
+				<Torsion k="10.46" p1="15" p2="14" p3="16" p4="17" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="8.368" p1="15" p2="14" p3="16" p4="17" periodicity="1" phase="0"/>
+				<Torsion k="10.46" p1="15" p2="14" p3="16" p4="18" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="43.932" p1="1" p2="6" p3="4" p4="5" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="4.6024" p1="4" p2="8" p3="6" p4="7" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="43.932" p1="8" p2="16" p3="14" p4="15" periodicity="2" phase="3.141592653589793"/>
+				<Torsion k="4.6024" p1="14" p2="18" p3="16" p4="17" periodicity="2" phase="3.141592653589793"/>
+			</Torsions>
+		</Force>
+		<Force alpha="0" cutoff="1" dispersionCorrection="1" ewaldTolerance=".0005" exceptionsUsePeriodic="0" forceGroup="0" includeDirectSpace="1" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" method="0" name="NonbondedForce" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="4">
+			<GlobalParameters/>
+			<ParticleOffsets/>
+			<ExceptionOffsets/>
+			<Particles>
+				<Particle eps=".06568879999999999" q=".1123" sig=".2649532787749369"/>
+				<Particle eps=".4577296" q="-.3662" sig=".3399669508423535"/>
+				<Particle eps=".06568879999999999" q=".1123" sig=".2649532787749369"/>
+				<Particle eps=".06568879999999999" q=".1123" sig=".2649532787749369"/>
+				<Particle eps=".359824" q=".5972" sig=".3399669508423535"/>
+				<Particle eps=".87864" q="-.5679" sig=".2959921901149463"/>
+				<Particle eps=".7112800000000001" q="-.4157" sig=".3249998523775958"/>
+				<Particle eps=".06568879999999999" q=".2719" sig=".10690784617684071"/>
+				<Particle eps=".4577296" q=".0337" sig=".3399669508423535"/>
+				<Particle eps=".06568879999999999" q=".0823" sig=".2471353044121301"/>
+				<Particle eps=".4577296" q="-.1825" sig=".3399669508423535"/>
+				<Particle eps=".06568879999999999" q=".0603" sig=".2649532787749369"/>
+				<Particle eps=".06568879999999999" q=".0603" sig=".2649532787749369"/>
+				<Particle eps=".06568879999999999" q=".0603" sig=".2649532787749369"/>
+				<Particle eps=".359824" q=".5973" sig=".3399669508423535"/>
+				<Particle eps=".87864" q="-.5679" sig=".2959921901149463"/>
+				<Particle eps=".7112800000000001" q="-.4157" sig=".3249998523775958"/>
+				<Particle eps=".06568879999999999" q=".2719" sig=".10690784617684071"/>
+				<Particle eps=".4577296" q="-.149" sig=".3399669508423535"/>
+				<Particle eps=".06568879999999999" q=".0976" sig=".2471353044121301"/>
+				<Particle eps=".06568879999999999" q=".0976" sig=".2471353044121301"/>
+				<Particle eps=".06568879999999999" q=".0976" sig=".2471353044121301"/>
+			</Particles>
+			<Exceptions>
+				<Exception eps="0" p1="0" p2="1" q="0" sig="1"/>
+				<Exception eps="0" p1="0" p2="2" q="0" sig="1"/>
+				<Exception eps="0" p1="1" p2="2" q="0" sig="1"/>
+				<Exception eps="0" p1="0" p2="3" q="0" sig="1"/>
+				<Exception eps="0" p1="1" p2="3" q="0" sig="1"/>
+				<Exception eps="0" p1="2" p2="3" q="0" sig="1"/>
+				<Exception eps="0" p1="0" p2="4" q="0" sig="1"/>
+				<Exception eps="0" p1="1" p2="4" q="0" sig="1"/>
+				<Exception eps="0" p1="2" p2="4" q="0" sig="1"/>
+				<Exception eps="0" p1="3" p2="4" q="0" sig="1"/>
+				<Exception eps=".12012161257658839" p1="0" p2="5" q="-.053145975" sig=".2804727344449416"/>
+				<Exception eps="0" p1="1" p2="5" q="0" sig="1"/>
+				<Exception eps=".12012161257658839" p1="2" p2="5" q="-.053145975" sig=".2804727344449416"/>
+				<Exception eps=".12012161257658839" p1="3" p2="5" q="-.053145975" sig=".2804727344449416"/>
+				<Exception eps="0" p1="4" p2="5" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="0" p2="6" q="-.03890259166666667" sig=".2949765655762664"/>
+				<Exception eps="0" p1="1" p2="6" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="2" p2="6" q="-.03890259166666667" sig=".2949765655762664"/>
+				<Exception eps=".10807766844265286" p1="3" p2="6" q="-.03890259166666667" sig=".2949765655762664"/>
+				<Exception eps="0" p1="4" p2="6" q="0" sig="1"/>
+				<Exception eps="0" p1="5" p2="6" q="0" sig="1"/>
+				<Exception eps=".08670021359327784" p1="1" p2="7" q="-.08297481666666667" sig=".2234373985095971"/>
+				<Exception eps="0" p1="4" p2="7" q="0" sig="1"/>
+				<Exception eps=".12012161257658839" p1="5" p2="7" q="-.128676675" sig=".20145001814589353"/>
+				<Exception eps="0" p1="6" p2="7" q="0" sig="1"/>
+				<Exception eps=".2288648" p1="1" p2="8" q="-.010284116666666667" sig=".3399669508423535"/>
+				<Exception eps="0" p1="4" p2="8" q="0" sig="1"/>
+				<Exception eps=".3170881327580709" p1="5" p2="8" q="-.015948525" sig=".3179795704786499"/>
+				<Exception eps="0" p1="6" p2="8" q="0" sig="1"/>
+				<Exception eps="0" p1="7" p2="8" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="4" p2="9" q=".040957966666666665" sig=".2935511276272418"/>
+				<Exception eps="0" p1="6" p2="9" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="7" p2="9" q=".01864780833333333" sig=".17702157529448542"/>
+				<Exception eps="0" p1="8" p2="9" q="0" sig="1"/>
+				<Exception eps=".20291752979375635" p1="4" p2="10" q="-.09082416666666666" sig=".3399669508423535"/>
+				<Exception eps="0" p1="6" p2="10" q="0" sig="1"/>
+				<Exception eps=".08670021359327784" p1="7" p2="10" q="-.041351458333333334" sig=".2234373985095971"/>
+				<Exception eps="0" p1="8" p2="10" q="0" sig="1"/>
+				<Exception eps="0" p1="9" p2="10" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="6" p2="11" q="-.020888925000000003" sig=".2949765655762664"/>
+				<Exception eps="0" p1="8" p2="11" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="9" p2="11" q=".004135575" sig=".2560442915935335"/>
+				<Exception eps="0" p1="10" p2="11" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="6" p2="12" q="-.020888925000000003" sig=".2949765655762664"/>
+				<Exception eps="0" p1="8" p2="12" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="9" p2="12" q=".004135575" sig=".2560442915935335"/>
+				<Exception eps="0" p1="10" p2="12" q="0" sig="1"/>
+				<Exception eps="0" p1="11" p2="12" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="6" p2="13" q="-.020888925000000003" sig=".2949765655762664"/>
+				<Exception eps="0" p1="8" p2="13" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="9" p2="13" q=".004135575" sig=".2560442915935335"/>
+				<Exception eps="0" p1="10" p2="13" q="0" sig="1"/>
+				<Exception eps="0" p1="11" p2="13" q="0" sig="1"/>
+				<Exception eps="0" p1="12" p2="13" q="0" sig="1"/>
+				<Exception eps=".179912" p1="4" p2="14" q=".29725630000000003" sig=".3399669508423535"/>
+				<Exception eps="0" p1="6" p2="14" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="7" p2="14" q=".135338225" sig=".2234373985095971"/>
+				<Exception eps="0" p1="8" p2="14" q="0" sig="1"/>
+				<Exception eps="0" p1="9" p2="14" q="0" sig="1"/>
+				<Exception eps="0" p1="10" p2="14" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="11" p2="14" q=".030014325000000005" sig=".3024601148086452"/>
+				<Exception eps=".07687068162049819" p1="12" p2="14" q=".030014325000000005" sig=".3024601148086452"/>
+				<Exception eps=".07687068162049819" p1="13" p2="14" q=".030014325000000005" sig=".3024601148086452"/>
+				<Exception eps=".3952717606912996" p1="6" p2="15" q=".196730025" sig=".31049602124627107"/>
+				<Exception eps="0" p1="8" p2="15" q="0" sig="1"/>
+				<Exception eps=".12012161257658839" p1="9" p2="15" q="-.038948474999999996" sig=".2715637472635382"/>
+				<Exception eps=".3170881327580709" p1="10" p2="15" q=".086368125" sig=".3179795704786499"/>
+				<Exception eps="0" p1="14" p2="15" q="0" sig="1"/>
+				<Exception eps=".35564000000000007" p1="6" p2="16" q=".14400540833333336" sig=".3249998523775958"/>
+				<Exception eps="0" p1="8" p2="16" q="0" sig="1"/>
+				<Exception eps=".10807766844265286" p1="9" p2="16" q="-.028510091666666668" sig=".28606757839486296"/>
+				<Exception eps=".28529542140034425" p1="10" p2="16" q=".06322104166666667" sig=".3324834016099747"/>
+				<Exception eps="0" p1="14" p2="16" q="0" sig="1"/>
+				<Exception eps="0" p1="15" p2="16" q="0" sig="1"/>
+				<Exception eps=".08670021359327784" p1="8" p2="17" q=".007635858333333333" sig=".2234373985095971"/>
+				<Exception eps="0" p1="14" p2="17" q="0" sig="1"/>
+				<Exception eps=".12012161257658839" p1="15" p2="17" q="-.128676675" sig=".20145001814589353"/>
+				<Exception eps="0" p1="16" p2="17" q="0" sig="1"/>
+				<Exception eps=".2288648" p1="8" p2="18" q="-.004184416666666667" sig=".3399669508423535"/>
+				<Exception eps="0" p1="14" p2="18" q="0" sig="1"/>
+				<Exception eps=".3170881327580709" p1="15" p2="18" q=".07051425" sig=".3179795704786499"/>
+				<Exception eps="0" p1="16" p2="18" q="0" sig="1"/>
+				<Exception eps="0" p1="17" p2="18" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="14" p2="19" q=".04858040000000001" sig=".2935511276272418"/>
+				<Exception eps="0" p1="16" p2="19" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="17" p2="19" q=".022114533333333335" sig=".17702157529448542"/>
+				<Exception eps="0" p1="18" p2="19" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="14" p2="20" q=".04858040000000001" sig=".2935511276272418"/>
+				<Exception eps="0" p1="16" p2="20" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="17" p2="20" q=".022114533333333335" sig=".17702157529448542"/>
+				<Exception eps="0" p1="18" p2="20" q="0" sig="1"/>
+				<Exception eps="0" p1="19" p2="20" q="0" sig="1"/>
+				<Exception eps=".07687068162049819" p1="14" p2="21" q=".04858040000000001" sig=".2935511276272418"/>
+				<Exception eps="0" p1="16" p2="21" q="0" sig="1"/>
+				<Exception eps=".032844399999999996" p1="17" p2="21" q=".022114533333333335" sig=".17702157529448542"/>
+				<Exception eps="0" p1="18" p2="21" q="0" sig="1"/>
+				<Exception eps="0" p1="19" p2="21" q="0" sig="1"/>
+				<Exception eps="0" p1="20" p2="21" q="0" sig="1"/>
+			</Exceptions>
+		</Force>
+		<Force forceGroup="0" frequency="1" name="CMMotionRemover" type="CMMotionRemover" version="1"/>
+		<Force forceGroup="0" name="HarmonicAngleForce" type="HarmonicAngleForce" usesPeriodic="0" version="2">
+			<Angles>
+				<Angle a="1.911135530933791" k="292.88" p1="0" p2="1" p3="2"/>
+				<Angle a="1.911135530933791" k="292.88" p1="0" p2="1" p3="3"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="0" p2="1" p3="4"/>
+				<Angle a="2.101376419401173" k="669.44" p1="1" p2="4" p3="5"/>
+				<Angle a="2.035053907825388" k="585.76" p1="1" p2="4" p3="6"/>
+				<Angle a="1.911135530933791" k="292.88" p1="2" p2="1" p3="3"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="2" p2="1" p3="4"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="3" p2="1" p3="4"/>
+				<Angle a="2.0943951023931953" k="418.40000000000003" p1="4" p2="6" p3="7"/>
+				<Angle a="2.1275563581810877" k="418.40000000000003" p1="4" p2="6" p3="8"/>
+				<Angle a="2.1450096507010312" k="669.44" p1="5" p2="4" p3="6"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="6" p2="8" p3="9"/>
+				<Angle a="1.9146261894377796" k="669.44" p1="6" p2="8" p3="10"/>
+				<Angle a="1.9216075064457567" k="527.184" p1="6" p2="8" p3="14"/>
+				<Angle a="2.0601866490541068" k="418.40000000000003" p1="7" p2="6" p3="8"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="8" p2="10" p3="11"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="8" p2="10" p3="12"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="8" p2="10" p3="13"/>
+				<Angle a="2.101376419401173" k="669.44" p1="8" p2="14" p3="15"/>
+				<Angle a="2.035053907825388" k="585.76" p1="8" p2="14" p3="16"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="9" p2="8" p3="10"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="9" p2="8" p3="14"/>
+				<Angle a="1.9390607989657" k="527.184" p1="10" p2="8" p3="14"/>
+				<Angle a="1.911135530933791" k="292.88" p1="11" p2="10" p3="12"/>
+				<Angle a="1.911135530933791" k="292.88" p1="11" p2="10" p3="13"/>
+				<Angle a="1.911135530933791" k="292.88" p1="12" p2="10" p3="13"/>
+				<Angle a="2.0943951023931953" k="418.40000000000003" p1="14" p2="16" p3="17"/>
+				<Angle a="2.1275563581810877" k="418.40000000000003" p1="14" p2="16" p3="18"/>
+				<Angle a="2.1450096507010312" k="669.44" p1="15" p2="14" p3="16"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="16" p2="18" p3="19"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="16" p2="18" p3="20"/>
+				<Angle a="1.911135530933791" k="418.40000000000003" p1="16" p2="18" p3="21"/>
+				<Angle a="2.0601866490541068" k="418.40000000000003" p1="17" p2="16" p3="18"/>
+				<Angle a="1.911135530933791" k="292.88" p1="19" p2="18" p3="20"/>
+				<Angle a="1.911135530933791" k="292.88" p1="19" p2="18" p3="21"/>
+				<Angle a="1.911135530933791" k="292.88" p1="20" p2="18" p3="21"/>
+			</Angles>
+		</Force>
+	</Forces>
+</System>
diff --git a/examples/template.yaml b/examples/template.yaml
index baed0ca..1501906 100644
--- a/examples/template.yaml
+++ b/examples/template.yaml
@@ -21,7 +21,7 @@ system:
 # ---- Engine: the MD backend and thermodynamic settings ---------------------
 engine:
   md_engine: openmm             # openmm | gromacs | amber
-  platform_name: CUDA           # CUDA | CPU | OpenCL | Reference  (OpenMM)
+  platform_name: CUDA           # use CPU if OpenMM has no registered CUDA platform
   precision: mixed              # mixed | single | double
   temperature: 300.0            # Kelvin
   pressure: 1.0                 # bar
@@ -56,7 +56,7 @@ spawning:
 # ---- CV space: fixed physical CVs or a learned latent space ----------------
 # space_mode: fixed | pca | tica | tvae | vampnet | spib | deep-tica | deep-lda
 space_mode: vampnet
-adaptive_feature_type: distances        # distances | fitted_coords | phi_psi
+adaptive_feature_type: distances        # distances | fitted_coords | phi_psi (AIB9-only)
 retrain_freq: 5                         # retrain cadence for retrain_policy: fixed
 retrain_policy: fixed                   # fixed | vamp_adaptive (retrain on VAMP-2 drop)
 # vamp_retrain_tol: 0.1                 # relative VAMP-2 drop that triggers a retrain
@@ -125,6 +125,7 @@ binning:
 # ---- Execution: where walkers run ------------------------------------------
 execution:
   backend: local                        # local | slurm | pbs
+  # walker_timeout: 3600                 # local: kill a walker after N seconds (hang guard)
   # --- scheduler settings (slurm/pbs) ---
   # partition: gpu
   # account: my_alloc
diff --git a/mkdocs.yml b/mkdocs.yml
index edfc451..b4233e7 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -22,9 +22,23 @@ nav:
   - Adaptive binning: binning.md
   - Analysis & plotting: analysis.md
   - Execution (workstation & HPC): execution.md
+  - CLI reference: cli.md
+  - API reference: api.md
+  - References: references.md
   - Tutorial — adaptive MSM: tutorials/adaptive_msm.md
   - Tutorial — notebook: notebook_tutorial.md
 
+plugins:
+  - search
+  - mkdocstrings:
+      handlers:
+        python:
+          options:
+            show_root_heading: true
+            show_source: false
+            members_order: source
+            docstring_style: numpy
+
 markdown_extensions:
   - admonition
   - pymdownx.superfences
diff --git a/pyproject.toml b/pyproject.toml
index 508aa2b..061af94 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,11 +4,20 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "autosampler"
-version = "2.0.0"
+dynamic = ["version"]
 description = "Adaptive collective-variable sampling framework for molecular dynamics."
+readme = "README.md"
+license = { file = "LICENSE" }
 requires-python = ">=3.10"
 authors = [
-  { name = "AutoSampler contributors" }
+  # TODO(release): confirm the full author list, affiliations, and ORCIDs for the
+  # publication (mirror these in CITATION.cff).
+  { name = "Dibyendu Maity", email = "dibyendumaity1999@bose.res.in" },
+  { name = "Suman Chakrabarty", email = "chakrabarty.suman@gmail.com" },
+]
+keywords = [
+  "molecular-dynamics", "adaptive-sampling", "markov-state-model",
+  "collective-variables", "enhanced-sampling", "openmm",
 ]
 dependencies = [
   "numpy>=1.23",
@@ -18,12 +27,11 @@ dependencies = [
   "pyyaml>=6.0",
   "scikit-learn>=1.2",
   "MDAnalysis>=2.5",
-  "openmm>=8.0",
   "torch>=2.0",
   "deeptime>=0.4",
 ]
 classifiers = [
-  "Development Status :: 3 - Alpha",
+  "Development Status :: 4 - Beta",
   "Intended Audience :: Science/Research",
   "License :: OSI Approved :: MIT License",
   "Operating System :: OS Independent",
@@ -35,7 +43,21 @@ classifiers = [
   "Topic :: Scientific/Engineering :: Bio-Informatics",
 ]
 
+[project.urls]
+Homepage = "https://github.com/TeamSuman/AutoSampler"
+Repository = "https://github.com/TeamSuman/AutoSampler"
+Documentation = "https://github.com/TeamSuman/AutoSampler/tree/main/docs"
+Issues = "https://github.com/TeamSuman/AutoSampler/issues"
+Changelog = "https://github.com/TeamSuman/AutoSampler/blob/main/CHANGELOG.md"
+
 [project.optional-dependencies]
+# The default OpenMM engine. Kept optional because OpenMM is most reliably
+# installed from conda-forge; the engine is imported lazily so the base package
+# installs and imports without it. `pip install 'autosampler[openmm]'` pulls the
+# PyPI wheel where available.
+openmm = [
+  "openmm>=8.0",
+]
 deep-tica = [
   "lightning>=2.0",
   "mlcolvar>=1.0",
@@ -51,7 +73,7 @@ docs = [
   "mkdocs-material>=9.0",
 ]
 all = [
-  "autosampler[deep-tica,examples,test,docs]",
+  "autosampler[openmm,deep-tica,examples,test,docs]",
 ]
 
 [tool.setuptools.packages.find]
@@ -62,6 +84,9 @@ exclude = ["examples*", "tests*", "latex*"]
 [tool.setuptools]
 include-package-data = false
 
+[tool.setuptools.dynamic]
+version = { attr = "autosampler.__version__" }  # single source: autosampler/__init__.py
+
 [project.scripts]
 autosampler = "autosampler.cli:main"
 autosampler-run = "autosampler.cli:main"
@@ -97,3 +122,5 @@ ignore = [
 [tool.ruff.lint.per-file-ignores]
 "__init__.py" = ["F401"]  # re-exports
 "tests/*" = ["E402"]  # importorskip before imports
+"autosampler/core.py" = ["E402"]  # warnings.filterwarnings must precede imports
+"autosampler/engines/openmm.py" = ["F403", "F405"]  # `from openmm import *` API
diff --git a/tests/test_checkpoint_delta.py b/tests/test_checkpoint_delta.py
new file mode 100644
index 0000000..8171253
--- /dev/null
+++ b/tests/test_checkpoint_delta.py
@@ -0,0 +1,79 @@
+"""Regression tests for delta-checkpointed history (save/reconstruct/resume).
+
+History is written incrementally: each ``iter_*/history.pkl`` holds only the
+entries since the previous checkpoint, and ``load`` / ``paths.load_history``
+reconstruct the full history by merging the deltas. These tests cover the
+round-trip, checkpoint-frequency gaps, corruption tolerance, and atomicity.
+"""
+
+from __future__ import annotations
+
+import pickle
+
+import pytest
+
+from autosampler.checkpoints.manager import CheckpointManager
+
+
+def test_delta_history_roundtrip(tmp_path):
+    mgr = CheckpointManager(str(tmp_path))
+    history: dict = {}
+    for it in range(4):
+        history[it] = {"frames": [{"iter": it}]}
+        mgr.save(it, None, {"s": it}, {"b": it}, dict(history))
+
+    # Each delta holds exactly the one new key for that iteration.
+    for it in range(4):
+        with open(tmp_path / f"iter_{it}" / "history.pkl", "rb") as f:
+            assert set(pickle.load(f)) == {it}
+
+    _, scaler, bin_state, full, _ = mgr.load(3)
+    assert set(full) == {0, 1, 2, 3}
+    assert scaler == {"s": 3} and bin_state == {"b": 3}
+
+
+def test_delta_history_with_checkpoint_gaps(tmp_path):
+    # checkpoint_freq > 1: only iters 0 and 2 are written, but history accrued
+    # every iteration must be fully reconstructed.
+    mgr = CheckpointManager(str(tmp_path))
+    mgr.save(0, None, {}, {}, {0: "a"})
+    mgr.save(2, None, {}, {}, {0: "a", 1: "b", 2: "c"})
+
+    with open(tmp_path / "iter_2" / "history.pkl", "rb") as f:
+        assert set(pickle.load(f)) == {1, 2}  # delta covers the gap
+
+    _, _, _, full, _ = mgr.load(2)
+    assert full == {0: "a", 1: "b", 2: "c"}
+
+
+def test_delta_history_tolerates_corrupt_delta(tmp_path):
+    mgr = CheckpointManager(str(tmp_path))
+    mgr.save(0, None, {}, {}, {0: "a"})
+    mgr.save(1, None, {}, {}, {0: "a", 1: "b"})
+    # Simulate a crash-truncated earlier delta.
+    (tmp_path / "iter_0" / "history.pkl").write_bytes(b"\x80\x04truncated")
+
+    _, _, _, full, _ = mgr.load(1)  # must not raise
+    assert 1 in full  # the newer, intact delta is still recovered
+
+
+def test_checkpoint_writes_are_atomic(tmp_path):
+    mgr = CheckpointManager(str(tmp_path))
+    mgr.save(0, None, {"s": 0}, {"b": 0}, {0: "a"})
+    # No leftover temp files from the atomic write+replace.
+    assert not list(tmp_path.glob("iter_*/*.tmp"))
+
+
+def test_paths_load_history_reconstructs_full_history(tmp_path):
+    from autosampler.paths import load_history
+
+    mgr = CheckpointManager(str(tmp_path / "checkpoints"))
+    mgr.save(0, None, {}, {}, {0: {"frames": []}})
+    mgr.save(1, None, {}, {}, {0: {"frames": []}, 1: {"frames": []}})
+
+    # Reads the whole run dir; must merge deltas, not return the last window only.
+    assert set(load_history(tmp_path)) == {0, 1}
+    assert set(load_history(tmp_path, checkpoint=1)) == {0, 1}
+
+    with pytest.raises(FileNotFoundError):
+        load_history(tmp_path, checkpoint=99)
diff --git a/tests/test_execution.py b/tests/test_execution.py
index a2b681a..03e5f21 100644
--- a/tests/test_execution.py
+++ b/tests/test_execution.py
@@ -109,6 +109,36 @@ def test_local_backend_runs_all_walkers(tmp_path, monkeypatch):
     assert results == [True, False, True]
 
 
+class _HangingExecutor:
+    """Hands out futures that never complete, to simulate a hung walker."""
+
+    def __init__(self, *a, **k):
+        self._processes = {}
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *a):
+        return False
+
+    def submit(self, fn, *args):
+        return Future()  # never set → never done
+
+
+def test_local_backend_walker_timeout(tmp_path, monkeypatch):
+    import time
+
+    monkeypatch.setattr(local_mod, "ProcessPoolExecutor", _HangingExecutor)
+    backend = ExecutionBackendFactory.get(
+        "local", gpu_ids=[0], max_workers=2, walker_timeout=0.3
+    )
+    start = time.monotonic()
+    results = backend.execute(_tasks(tmp_path, ["hang", "hang"]))
+    elapsed = time.monotonic() - start
+    assert results == [False, False]  # timed-out batch reported as failed
+    assert elapsed < 5.0  # and it did not hang
+
+
 # ── fake scheduler command runner ───────────────────────────────────────────
 def _fake_runner(submit_id="4242"):
     def runner(cmd, timeout):
diff --git a/tests/test_reproducibility.py b/tests/test_reproducibility.py
new file mode 100644
index 0000000..bd4f79e
--- /dev/null
+++ b/tests/test_reproducibility.py
@@ -0,0 +1,56 @@
+"""Reproducibility: seed plumbing and deterministic CV training."""
+
+from __future__ import annotations
+
+import warnings
+
+import numpy as np
+import pytest
+
+warnings.filterwarnings("ignore")
+
+
+def test_seed_manager_enables_deterministic_mode():
+    from autosampler.utils.seeds import SeedManager
+
+    SeedManager(123).set_seed()  # must not raise even with deterministic algos on
+    a = np.random.rand(5)
+    SeedManager(123).set_seed()
+    b = np.random.rand(5)
+    np.testing.assert_array_equal(a, b)
+
+
+def test_adaptive_model_carries_seed():
+    torch = pytest.importorskip("torch")
+    pytest.importorskip("deeptime")
+    from autosampler.spaces.model import AdaptiveSpaceModel
+
+    m = AdaptiveSpaceModel(space_mode="vampnet", seed=7)
+    assert m.seed == 7
+
+    # Default seed is defined (round-trips through ensure_config_defaults).
+    m2 = AdaptiveSpaceModel(space_mode="vampnet")
+    m2.ensure_config_defaults()
+    assert isinstance(m2.seed, int)
+    del torch
+
+
+def test_same_seed_gives_identical_vampnet_projection():
+    pytest.importorskip("torch")
+    pytest.importorskip("deeptime")
+    from autosampler.spaces.model import AdaptiveSpaceModel
+
+    rng = np.random.default_rng(0)
+    # Two metastable blobs, ordered by walker then time (2 walkers).
+    feats = np.vstack(
+        [rng.normal(-1, 0.1, (50, 4)), rng.normal(1, 0.1, (50, 4))]
+    ).astype(np.float32)
+
+    def project():
+        m = AdaptiveSpaceModel(
+            space_mode="vampnet", seed=42, lagtime=2, latent_dim=2, epochs=3
+        )
+        m.fit(feats, walker_length=50, n_walkers=2)
+        return np.asarray(m.project(feats))
+
+    np.testing.assert_allclose(project(), project(), rtol=1e-5, atol=1e-5)
diff --git a/tests/test_spawner_weighted_choice.py b/tests/test_spawner_weighted_choice.py
new file mode 100644
index 0000000..5a50179
--- /dev/null
+++ b/tests/test_spawner_weighted_choice.py
@@ -0,0 +1,22 @@
+"""Regression: density/voronoi weighted choice must not crash when most weights
+are zero (target mode zeroes out distant bins)."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from autosampler.spawners.density import _weighted_choice
+
+
+def test_weighted_choice_with_mostly_zero_weights():
+    rows = np.arange(4)
+    weights = np.array([1.0, 0.0, 0.0, 0.0])  # only one nonzero
+    out = _weighted_choice(rows, weights, top_n=3)  # would raise pre-fix
+    assert len(out) == 3
+    assert set(out).issubset(set(rows.tolist()))
+
+
+def test_weighted_choice_all_zero_falls_back_to_uniform():
+    rows = np.arange(5)
+    out = _weighted_choice(rows, np.zeros(5), top_n=2)
+    assert len(out) == 2 and set(out).issubset(set(rows.tolist()))
diff --git a/tests/test_spawners_and_paths.py b/tests/test_spawners_and_paths.py
new file mode 100644
index 0000000..d78e6eb
--- /dev/null
+++ b/tests/test_spawners_and_paths.py
@@ -0,0 +1,110 @@
+"""Coverage for the previously-untested spawners (voronoi/lof/fps) and the
+trajectory-lineage helpers in autosampler.paths."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from autosampler.spawners.base import SpawnerFactory
+
+
+def _cloud(n=40, seed=0):
+    rng = np.random.default_rng(seed)
+    return np.vstack([rng.normal(-2, 0.3, (n // 2, 2)), rng.normal(2, 0.3, (n // 2, 2))])
+
+
+@pytest.mark.parametrize(
+    "scheme,kwargs",
+    [
+        ("fps", {}),
+        ("lof", {"n_neighbors": 8}),
+        ("voronoi", {"n_clusters": 10}),
+    ],
+)
+def test_spawner_explore_returns_valid_indices(scheme, kwargs):
+    if scheme == "voronoi":
+        pytest.importorskip("shapely")
+    if scheme == "lof":
+        pytest.importorskip("sklearn")
+    points = _cloud()
+    spawner = SpawnerFactory.get(scheme, mode="explore", **kwargs)
+    idx = spawner.sample(points, top_n=6)
+    assert len(idx) == 6
+    assert all(0 <= int(i) < len(points) for i in idx)
+
+
+@pytest.mark.parametrize("scheme", ["fps", "lof", "voronoi"])
+def test_spawner_target_mode(scheme):
+    if scheme == "voronoi":
+        pytest.importorskip("shapely")
+    if scheme == "lof":
+        pytest.importorskip("sklearn")
+    points = _cloud()
+    kwargs = {"target": [2.0, 2.0]}
+    if scheme == "lof":
+        kwargs["n_neighbors"] = 8
+    if scheme == "voronoi":
+        kwargs["n_clusters"] = 10
+    spawner = SpawnerFactory.get(scheme, mode="target", **kwargs)
+    idx = spawner.sample(points, top_n=5)
+    assert len(idx) == 5 and all(0 <= int(i) < len(points) for i in idx)
+
+
+def test_spawner_factory_unknown_scheme():
+    with pytest.raises((ValueError, KeyError)):
+        SpawnerFactory.get("nope")
+
+
+# ── paths / lineage helpers ─────────────────────────────────────────────────
+from autosampler import paths  # noqa: E402
+
+
+def test_frame_key_and_frameref_roundtrip():
+    assert paths.frame_key(3, 1, 4) == "3:1:4"
+    ref = paths.FrameRef(
+        iteration=1, walker=0, frame=2, trajectory="t.xtc", cv=(0.1, 0.2), parent=None
+    )
+    again = paths.FrameRef.from_dict(ref.to_dict())
+    assert again == ref
+    assert again.key == "1:0:2"
+
+
+def _records():
+    # 2 walkers x 2 frames, non-xtc so expected_frames is used (no real I/O).
+    pts = np.array([[0.0, 0.0], [1.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
+    return paths.build_frame_records(
+        iteration=0,
+        trajectories=["w0.dcd", "w1.dcd"],
+        points=pts,
+        walker_parents=[None, None],
+        expected_frames=2,
+    )
+
+
+def test_build_frame_records_and_lineage():
+    recs = _records()
+    assert len(recs) == 4
+    # within-walker parent links: frame 1 points at frame 0 of the same walker.
+    assert recs[1]["parent"] == paths.frame_key(0, 0, 0)
+    assert recs[0]["parent"] is None
+
+    mapped = paths.map_global_frame(recs, 2)
+    assert mapped["walker"] == 1 and mapped["frame"] == 0
+    with pytest.raises(IndexError):
+        paths.map_global_frame(recs, 99)
+
+
+def test_history_records_and_nearest():
+    history = {0: {"frames": _records()}}
+    refs = paths.history_records(history)
+    assert len(refs) == 4
+    near = paths.nearest_record(refs, np.array([2.9, 2.9]))
+    assert near.cv == (3.0, 3.0)
+    with pytest.raises(ValueError):
+        paths.nearest_record(refs, np.array([0.0, 0.0, 0.0]))  # dim mismatch
+
+
+def test_history_records_empty_raises():
+    with pytest.raises(ValueError):
+        paths.history_records({0: {"frames": []}})
diff --git a/tools/generate_review_pdf.py b/tools/generate_review_pdf.py
new file mode 100644
index 0000000..b21240a
--- /dev/null
+++ b/tools/generate_review_pdf.py
@@ -0,0 +1,352 @@
+"""Generate a structured PDF of the AutoSampler production-readiness review."""
+
+import os
+
+from reportlab.lib import colors
+from reportlab.lib.enums import TA_LEFT
+from reportlab.lib.pagesizes import A4
+from reportlab.lib.styles import ParagraphStyle, getSampleStyleSheet
+from reportlab.lib.units import mm
+from reportlab.platypus import (
+    HRFlowable,
+    ListFlowable,
+    ListItem,
+    PageBreak,
+    Paragraph,
+    SimpleDocTemplate,
+    Spacer,
+    Table,
+    TableStyle,
+)
+
+OUT = os.path.join(
+    os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
+    "AutoSampler_Review.pdf",
+)
+
+NAVY = colors.HexColor("#1A237E")
+INDIGO = colors.HexColor("#3949AB")
+RED = colors.HexColor("#B71C1C")
+AMBER = colors.HexColor("#E65100")
+GREEN = colors.HexColor("#2E7D32")
+LGREY = colors.HexColor("#ECEFF1")
+DGREY = colors.HexColor("#37474F")
+
+ss = getSampleStyleSheet()
+H1 = ParagraphStyle("H1", parent=ss["Heading1"], textColor=NAVY, fontSize=15,
+                    spaceBefore=12, spaceAfter=6)
+H2 = ParagraphStyle("H2", parent=ss["Heading2"], textColor=INDIGO, fontSize=12,
+                    spaceBefore=9, spaceAfter=3)
+BODY = ParagraphStyle("Body", parent=ss["BodyText"], fontSize=9.5, leading=13.5,
+                      alignment=TA_LEFT, spaceAfter=4)
+SMALL = ParagraphStyle("Small", parent=BODY, fontSize=8.4, textColor=DGREY)
+TITLE = ParagraphStyle("Title", parent=ss["Title"], textColor=NAVY, fontSize=25,
+                       leading=29, spaceAfter=6)
+SUB = ParagraphStyle("Sub", parent=ss["Title"], textColor=INDIGO, fontSize=13,
+                     leading=17, spaceAfter=4)
+CELL = ParagraphStyle("Cell", parent=BODY, fontSize=8.3, leading=10.8, spaceAfter=0)
+CELLH = ParagraphStyle("CellH", parent=CELL, textColor=colors.white,
+                       fontName="Helvetica-Bold")
+
+story = []
+
+
+def bullets(items, style=BODY):
+    return ListFlowable(
+        [ListItem(Paragraph(t, style), leftIndent=10, value="•") for t in items],
+        bulletType="bullet", start="•", leftIndent=12, spaceBefore=1, spaceAfter=6,
+    )
+
+
+def rule():
+    story.append(Spacer(1, 3))
+    story.append(HRFlowable(width="100%", thickness=0.6, color=INDIGO))
+    story.append(Spacer(1, 5))
+
+
+def sev_table(rows, col_widths, header_bg=NAVY):
+    data = [[Paragraph(c, CELLH if i == 0 else CELL) for c in r]
+            for i, r in enumerate(rows)]
+    tbl = Table(data, colWidths=col_widths, repeatRows=1)
+    style = [
+        ("BACKGROUND", (0, 0), (-1, 0), header_bg),
+        ("ROWBACKGROUNDS", (0, 1), (-1, -1), [colors.white, LGREY]),
+        ("GRID", (0, 0), (-1, -1), 0.4, colors.HexColor("#B0BEC5")),
+        ("VALIGN", (0, 0), (-1, -1), "TOP"),
+        ("LEFTPADDING", (0, 0), (-1, -1), 5), ("RIGHTPADDING", (0, 0), (-1, -1), 5),
+        ("TOPPADDING", (0, 0), (-1, -1), 4), ("BOTTOMPADDING", (0, 0), (-1, -1), 4),
+    ]
+    tbl.setStyle(TableStyle(style))
+    story.append(tbl)
+
+
+# ---------------- Title page ----------------
+story.append(Spacer(1, 36))
+story.append(Paragraph("AutoSampler", TITLE))
+story.append(Paragraph("Production &amp; Publication Readiness Review", SUB))
+rule()
+story.append(Paragraph(
+    "An independent code review of the <b>main</b> branch, conducted from the dual "
+    "perspective of a code reviewer and a prospective user, ahead of a research "
+    "publication announcing the package.", BODY))
+story.append(Spacer(1, 4))
+story.append(Paragraph(
+    "Review date: 2026-06-30 &nbsp;|&nbsp; Method: five parallel domain reviewers "
+    "(code correctness, documentation, tests/CI/packaging, examples/UX, test-coverage "
+    "map) plus direct lead verification of all headline findings.", SMALL))
+story.append(Spacer(1, 10))
+
+story.append(Paragraph("Verdict", H2))
+story.append(Paragraph(
+    "<b>Strong scientific core; not yet release-ready as found.</b> The algorithmic "
+    "layer (MSM estimation &amp; convergence, adaptive binning, weighted ensemble, "
+    "VAMP-2 feature selection, retraining, scheduler logic) is well-engineered and "
+    "genuinely well-tested, and the docs site, annotated template, and rendered "
+    "notebook are excellent for a research package. However, the review found "
+    "<b>~6 release-gating blockers</b> and a band of correctness, reproducibility, "
+    "packaging, and example/test gaps. A subset of blockers has already been fixed "
+    "during this review (Section 2); the remainder is itemised with a prioritised "
+    "plan.", BODY))
+story.append(Spacer(1, 6))
+
+sev_table([
+    ["Area", "As-found state", "Severity band"],
+    ["Release gates", "Test suite RED; documented install broken; Quick Start example "
+     "could not run", "Blocker"],
+    ["Correctness / robustness", "Walker-failure &amp; hang handling, engine/scheduler "
+     "edge cases, target-mode crash", "Major"],
+    ["Reproducibility", "Seeding not actually deterministic across the training loop "
+     "(matters for a paper)", "Major"],
+    ["Packaging / release", "Incomplete PyPI metadata, heavy deps, partial CI lint, no "
+     "release/DOI", "Major"],
+    ["Docs / examples / tests", "Gaps: API ref, citations, several features undemoed, "
+     "orchestration untested", "Major / Minor"],
+], [34 * mm, 104 * mm, 32 * mm])
+story.append(Spacer(1, 8))
+story.append(Paragraph(
+    "Scope note: findings are rated Blocker / Major / Minor / Nit. CONFIRMED items "
+    "were traced in source by the lead reviewer; others are reviewer-reported with a "
+    "code citation. One widely-suspected bug was investigated and dismissed as a "
+    "false positive (see Section 6).", SMALL))
+story.append(PageBreak())
+
+# ---------------- Section 1: strengths ----------------
+story.append(Paragraph("1. What is already strong", H1))
+rule()
+story.append(bullets([
+    "<b>MSM subsystem</b> (estimator + convergence monitor): connected-set "
+    "restriction, MLE/Bayesian estimation, implied timescales, VAMP-2, PCCA+, and a "
+    "flux-weighted transition-matrix convergence criterion — with real, "
+    "behaviour-asserting tests (recovers 3-state systems, serialization round-trips).",
+    "<b>Adaptive binning</b> (gradient / mab / eigenvector), <b>weighted ensemble</b> "
+    "(weight-conserving split/merge), <b>VAMP-2 feature selection</b>, and "
+    "<b>retraining policies</b> are cleanly factored and tested.",
+    "<b>Execution backends</b> (local / SLURM / PBS) behind a clean strategy "
+    "interface; scheduler logic is unit-tested with a fake command runner.",
+    "<b>Developer experience</b>: annotated input-file template + autosampler-init, a "
+    "full MkDocs site, and a rendered Jupyter notebook tutorial that runs on synthetic "
+    "data without a GPU.",
+    "Test quality where present is high — assertions check real numerical/behavioural "
+    "outcomes, not merely the absence of exceptions.",
+]))
+
+# ---------------- Section 2: fixed during review ----------------
+story.append(Paragraph("2. Fixed during this review", H1))
+rule()
+story.append(Paragraph(
+    "The following were implemented and pushed to branch "
+    "<font face='Courier'>claude/autosampler-analysis-plan-jj0n1q</font> (suite green "
+    "at 101 tests). They are recorded here so the report doubles as a change log.", BODY))
+sev_table([
+    ["#", "Item", "Resolution"],
+    ["1", "Test suite RED — examples/template.yaml drifted from the module template",
+     "Regenerated from the single source; suite green."],
+    ["2", "Documented conda install broken (env.yml pinned pydantic 1.x; missing "
+     "shapely)", "Bumped to pydantic&gt;=2.0, added shapely + deep-CV backends."],
+    ["3", "Delta-checkpoint resume truncated history (broke autosampler-path); "
+     "non-atomic writes", "Reconstruct full history across deltas; atomic writes; "
+     "tolerate a corrupt delta; added regression tests."],
+    ["4", "Local backend aborted the whole iteration if one walker failed",
+     "Catch per-walker failure, log, mark unsuccessful (matches scheduler path)."],
+    ["5", "Incomplete PyPI metadata; version duplicated",
+     "Added readme/license/urls/keywords/authors; single-sourced version; Beta status."],
+    ["6", "No citation infrastructure", "Added CITATION.cff + README 'How to cite' "
+     "(author list / ORCIDs / DOI flagged as TODO)."],
+    ["7", "AlaD Quick Start unrunnable (hardcoded dev path; silent platform typo)",
+     "Removed the path, fixed platform_name; documented the GROMACS FF requirement."],
+    ["8", "Doc drift (blob/devel link; phi_psi advertised as generic)",
+     "Fixed link to main; documented phi_psi as AIB9-specific."],
+], [8 * mm, 88 * mm, 74 * mm])
+story.append(PageBreak())
+
+# ---------------- Section 3: blockers ----------------
+story.append(Paragraph("3. Release-gating blockers", H1))
+rule()
+story.append(Paragraph(
+    "Must be resolved before a general release. Items 1–4, 6 below are addressed in "
+    "Section 2; the rest remain open.", BODY))
+sev_table([
+    ["Blocker", "Detail &amp; status"],
+    ["Red test suite",
+     "tests/test_input_template.py failed on main. <b>FIXED.</b>"],
+    ["Broken documented install",
+     "env.yml pydantic 1.x vs Pydantic-v2 code; missing shapely. <b>FIXED.</b>"],
+    ["Delta-checkpoint regression",
+     "autosampler-path read a truncated (delta-only) history; non-atomic writes. "
+     "<b>FIXED.</b>"],
+    ["Quick Start cannot run",
+     "AlaD needs an external GROMACS force field; no CPU-only hello-world exists. "
+     "<b>PARTIAL</b> — de-risked; a self-contained alanine-dipeptide example is the "
+     "agreed next step."],
+    ["phi_psi crashes off-AIB9",
+     "adaptive_feature_type: phi_psi hard-requires 9 AIB residues but is advertised "
+     "generically. <b>Documented;</b> code still AIB9-only — rename or generalise."],
+    ["No citation / DOI",
+     "Gates citability for the paper. <b>Scaffolded</b> (CITATION.cff); needs the "
+     "confirmed author list, ORCIDs, affiliations, and a Zenodo DOI."],
+], [34 * mm, 136 * mm])
+
+# ---------------- Section 4: major correctness ----------------
+story.append(Paragraph("4. Major — correctness &amp; robustness", H1))
+rule()
+sev_table([
+    ["Finding", "Risk", "Location"],
+    ["No MD timeout / watchdog on the default (OpenMM, in-process) engine or the local "
+     "backend; scheduler poll loop has no overall deadline", "A hung GPU/driver stalls "
+     "the campaign forever", "execution/local.py; engines/openmm.py; scheduler.py"],
+    ["OpenMM engine only returns True or raises; NaN-recovery re-steps outside any "
+     "guard", "An unstable spawn crashes the worker instead of reporting failure",
+     "engines/openmm.py"],
+    ["Subprocess engines move output after exit-0 without checking a trajectory was "
+     "produced", "mdrun exiting 0 with no/empty output → uncaught error",
+     "engines/gromacs.py; amber.py"],
+    ["Scheduler treats a failed poll command as 'job done'; SLURM job-id match is a "
+     "\\b-bounded substring (123 vs 1234)", "Transient squeue/qstat hiccup abandons a "
+     "healthy job; id collision", "execution/scheduler.py; slurm.py"],
+    ["expected_frames fallback + bare except in frame mapping; only the aggregate "
+     "count is checked", "Offsetting over/under-production mis-assigns CV rows → "
+     "corrupt lineage", "paths.py"],
+    ["Density/Voronoi _weighted_choice uses replace=False with many zero weights "
+     "(target mode)", "np.random.choice raises 'fewer non-zero entries than size'",
+     "spawners/density.py"],
+    ["Triclinic/truncated-octahedron boxes reduced to box diagonal + 90°", "Wrong PBC "
+     "for common Amber solvated systems, silently", "engines/amber.py; gromacs.py"],
+    ["deep-tica projection runs on CPU tensors while params may be on CUDA; no NaN/"
+     "convergence guard in CV training", "Device-mismatch RuntimeError; NaN CVs flow "
+     "silently into binning/spawning", "spaces/model.py; spib.py"],
+], [92 * mm, 44 * mm, 34 * mm])
+story.append(PageBreak())
+
+# ---------------- Section 5: reproducibility / packaging / docs ----------------
+story.append(Paragraph("5. Major — reproducibility, packaging, docs/examples/tests", H1))
+rule()
+
+story.append(Paragraph("Reproducibility (matters for the publication's claims)", H2))
+story.append(bullets([
+    "Seed is set once at startup, not before each CV retrain — 2nd+ retrains depend "
+    "on all intervening RNG draws.",
+    "torch.use_deterministic_algorithms is never enabled; tvae/vampnet DataLoaders get "
+    "no seeded generator (vampnet shuffles on the global RNG).",
+    "SPIB hardcodes seed=42 and resets the global torch RNG mid-loop, ignoring "
+    "random_seed and perturbing later draws.",
+    "Density/Voronoi/LOF/FPS spawners and the Voronoi binner draw from the unseeded "
+    "global np.random. Net: runs are not bit-reproducible despite the SeedManager.",
+]))
+
+story.append(Paragraph("Packaging / release (target: PyPI + conda)", H2))
+story.append(bullets([
+    "Heavy deps are mandatory; openmm is effectively conda-only, so a plain "
+    "pip install often fails to resolve — move openmm (and likely MDAnalysis) to "
+    "extras with a lazy import so the base install works.",
+    "CI lints only ~18 hand-picked files; the most-edited modules (core.py, config.py, "
+    "engines/*, checkpoints/manager.py) are unlinted (~250 ruff findings incl. "
+    "<font face='Courier'>from openmm import *</font> star imports). Lint the whole tree.",
+    "No coverage gate, no <font face='Courier'>mkdocs build --strict</font> check, no "
+    "tag-triggered release/Zenodo workflow; Python 3.12 advertised but untested.",
+    "Version single-sourcing &amp; full metadata: <b>DONE</b> in Section 2.",
+]))
+
+story.append(Paragraph("Docs, examples &amp; tests", H2))
+story.append(bullets([
+    "Config schema uses the Pydantic default (extra=ignore): a typo'd key passes "
+    "<font face='Courier'>--check</font> silently — consider extra=forbid (and fix the "
+    "example configs that rely on it).",
+    "configuration.md omits several real keys; CLI reference omits autosampler-init / "
+    "-analyze; no API reference (mkdocstrings); method citations are thin "
+    "(TICA/VAMPNet/PCCA+/MAB uncited).",
+    "No worked example for spib, deep-tica, lof, fps, we, target mode, pbs, or "
+    "mab/eigenvector binning — several are headline features.",
+    "core.py orchestration, MD engines, half the spawners (voronoi/lof/fps), paths.py, "
+    "and 3/5 CLIs are effectively untested; no end-to-end iteration test. (Delta "
+    "checkpointing is now covered.)",
+]))
+
+# ---------------- Section 6: minor + false positive ----------------
+story.append(Paragraph("6. Minor / nits &amp; a dismissed false positive", H1))
+rule()
+story.append(bullets([
+    "Leftover debug <font face='Courier'>print(\"EXACT SPAWN INDICES\", …)</font> "
+    "(core.py); hardcoded personal gmx path (paths.py); dead <font face='Courier'>hint"
+    "</font> variable (openmm.py).",
+    "FPS can return duplicate indices; TVAE BatchNorm + a final batch of size 1 raises; "
+    "binning.find_bins raises AttributeError if called before fit.",
+    "Broad <font face='Courier'>except Exception</font> swallowing without logging in "
+    "several spawner/binner fallbacks; unbounded scheduler _jobs/ artifact growth.",
+    "No NaN/inf validation of input CV points across spawners.",
+]))
+story.append(Spacer(1, 4))
+story.append(Paragraph(
+    "<b>Dismissed false positive.</b> A reviewer flagged the MSM least-counts weight "
+    "as <font face='Courier'>1/count²</font>. On tracing it, the per-frame value "
+    "correctly distributes the microstate-level <font face='Courier'>1/count</font> "
+    "least-counts weight across that microstate's frames — exactly mirroring the "
+    "MSM-guided path's <font face='Courier'>base / size_per_frame</font>. Left "
+    "unchanged.", SMALL))
+
+# ---------------- Section 7: plan ----------------
+story.append(Paragraph("7. Recommended path to release", H1))
+rule()
+story.append(Paragraph("Priority 1 — finish the blockers", H2))
+story.append(bullets([
+    "Ship the self-contained CPU-only alanine-dipeptide hello-world; generalise or "
+    "rename phi_psi; confirm the author list + mint a Zenodo DOI in CITATION.cff.",
+]))
+story.append(Paragraph("Priority 2 — correctness &amp; reproducibility", H2))
+story.append(bullets([
+    "Add an MD-walker timeout; make the OpenMM engine report failure instead of "
+    "raising; verify trajectories were produced; harden scheduler polling.",
+    "Plumb the configured seed through every fit / spawner / DataLoader and enable "
+    "torch determinism — required to claim reproducibility in the paper.",
+]))
+story.append(Paragraph("Priority 3 — packaging, docs, examples, tests", H2))
+story.append(bullets([
+    "Move openmm/MDAnalysis to extras with lazy imports; lint the whole tree in CI; "
+    "add coverage + a strict docs build + a release/Zenodo workflow; test on 3.12.",
+    "Add an API reference + a references page; complete the config &amp; CLI docs; add "
+    "examples for the undemoed features; broaden tests to orchestration, engines, the "
+    "remaining spawners, paths, and CLIs, plus one end-to-end iteration test.",
+]))
+story.append(Spacer(1, 8))
+story.append(Paragraph(
+    "Bottom line: the science and architecture are publication-grade; the gap to "
+    "release is engineering hardening, reproducibility, packaging, and coverage — all "
+    "tractable, and partly already underway on the review branch.", BODY))
+
+
+def footer(canvas, doc):
+    canvas.saveState()
+    canvas.setFont("Helvetica", 7.5)
+    canvas.setFillColor(DGREY)
+    canvas.drawString(20 * mm, 12 * mm, "AutoSampler — Production & Publication Readiness Review")
+    canvas.drawRightString(190 * mm, 12 * mm, f"Page {doc.page}")
+    canvas.restoreState()
+
+
+doc = SimpleDocTemplate(
+    OUT, pagesize=A4, leftMargin=18 * mm, rightMargin=18 * mm,
+    topMargin=16 * mm, bottomMargin=18 * mm,
+    title="AutoSampler Production & Publication Readiness Review",
+    author="AutoSampler review",
+)
+doc.build(story, onFirstPage=footer, onLaterPages=footer)
+print("WROTE", OUT)