Phase 10: conformer ensembles, transition structures in series, parallel channels

CHANGELOG.md

@@ -6,6 +6,72 @@ Notable changes to Kinisot. Format follows [Keep a Changelog](https://keepachang
 
 ### Added
 
+- Conformer ensembles (IMPLEMENTATION_PLAN.md, Phase 10). A species given
+  as several files, `kinisot --rct gs_1.out gs_2.out --ts ts_*.out`, or in
+  Python as a list or `Conformers(files, free_energies=None,
+  degeneracy=None)`, is an ensemble of conformers with the same atom
+  numbering, in fast equilibrium.
+  - KIE = Σ x_i ρ_i / Σ y_j ρ‡_j, with the weights of the light isotopologue,
+    which is exact; for a transition structure, the weight is its share of
+    the rate, including its tunnelling factor. EQEs work the same way.
+    docs/theory.md has a new section 7 with the derivation.
+  - `--weights` (`weights=`): `qrrho` (default; quasi-harmonic free
+    energies from Kinisot's frequencies with Grimme's entropy interpolation,
+    through GoodVibes' functions), `rrho`, `user`, `lowest` or `equal`.
+    `--energies TABLE` gives your own free energies and degeneracies.
+  - The result, `EnsembleIsotopeEffect`, lists each conformer's free
+    energy, share and own KIE, and gives the KIE of the lowest conformers,
+    the effective number of transition-structure conformers and the range
+    when each free energy moves by ±0.5 kcal/mol (`--weight-uncertainty`).
+    The text output prints the conformer table and `KIE (ensemble) @ T`
+    lines; JSON and CSV carry the same.
+  - Conformers of a species must have the same atoms in the same order
+    (an error names the first difference). Different bonding, duplicate
+    structures and different levels of theory are warnings.
+  - One conformer per species gives the ordinary `IsotopeEffect`, unchanged.
+  - `kinisot.equivalent_positions()` gives the exact isotope effect of
+    positions made equivalent by fast motion from one result per placement
+    of the label.
+  - Over the 18 Shi epoxidation transition structures, `compute_kie`
+    reproduces the prototype ensemble (`benchmarks/shi_epoxidation/ensemble`)
+    to 2 × 10⁻⁶ at every site. Benchmark cases may list conformers per
+    species.
+- Transition structures in series and parallel channels (Phase 10;
+  docs/theory.md, sections 7a and 7b).
+  - `compute_kie(rct, ts=Series([ts_1, ts_2, ...]), iso=...)` combines
+    steps none of which alone commits the substrate:
+    KIE = Σ w_n KIE_n with w_n ∝ exp(+G_n/RT)/κ_L,n. The weights come from a
+    commitment factor (two steps), the steps' free energies, or computed
+    free energies. The result gives each step's KIE and share, C_f, a
+    sensitivity range, `kie_at(C_f)` and `commitment_for(measured)`.
+  - `channels([dict(rct=..., ts=..., iso=...), ...])` combines parallel
+    routes with their own reactants, transition structures or labels:
+    1/KIE = Σ y_c/KIE_c. The shares are given (a measured selectivity),
+    built from barriers, or computed; `amounts` covers reactants that do not
+    interconvert, and `selectivity_for(measured)` inverts two channels.
+  - A step or a channel may be a conformer ensemble, and a channel may be a
+    series.
+  - `series_kie()` and `channel_kie()` apply the formulas to KIEs from
+    anywhere. They reproduce the Wittig predictions of Chen et al. (JACS
+    2014; 1.028 from the free energies, 1.033 from the trajectories) and
+    Figure 3d of van Dijk et al. (Nat. Catal. 2021; s = 3.3).
+  - Checks: the series formula matches the slowest eigenvalue of a
+    three-step rate matrix to 7 × 10⁻⁸, and channels computed from one
+    reactant reproduce the conformer ensemble of their transition
+    structures to 10⁻¹³.
+- Worked example `examples/conformers`: eight GFN2-xTB conformers of allyl
+  vinyl ether and the chair and boat Claisen transition structures, made by
+  `scripts/make_conformer_example.py` (conformer search by dihedral
+  rotation, duplicates and mirror images recognized, degeneracies from
+  chirality, each saddle point checked to connect reactant and product).
+  ¹³C and ¹⁷O KIEs hardly depend on the conformer; the C4-d₂ KIE runs from
+  0.926 to 0.984, and the ensemble (0.968) differs from the lowest pair by
+  0.010. `run_examples.sh` and `tests/test_examples.py` include it.
+- `kinisot --job job.json`: a JSON job file with the structures, labels and
+  settings, for series, channels, conformer ensembles with free energies,
+  and several isotopologues in one run (docs/job_files.md). The results
+  file prints the steps or channels with their shares, and `--json` and
+  `--csv` carry an `isotopologue` key.
 - Five benchmark cases from the SI of the PyQuiverHS paper (Grazioli, Ly,
   Sabetnejad, Mattapalli, Nguyen and O'Leary, ChemRxiv 2026), included with
   the authors' agreement. Each has its Gaussian files, the PyQuiverHS input
@@ -177,6 +243,14 @@ Notable changes to Kinisot. Format follows [Keep a Changelog](https://keepachang
 
 - The calculator registry records which keyword receives the `:model` part
   of a `--calc` specification (`method` for tblite, `model` elsewhere).
+- The benchmark runner averages equivalent positions (`iso_average`,
+  `reference_average`) exactly: the isotope ratios are averaged over the
+  placements of the label on each side, which for a difference on one side
+  is the harmonic mean of the separate KIEs. It took the geometric mean
+  before. Two reported values move: the Diels–Alder H3 KIE by 1 × 10⁻⁴
+  (0.9904 to 0.9905) and the singly labelled TS1B nitroarene KIE by
+  1 × 10⁻⁴ (1.0305 to 1.0304; 1.0321 to 1.0319 with Bell tunnelling). The
+  mean absolute deviation over the 34 measured positions stays 0.0058.
 
 ## [2.5.0] - 2026-09-25
 

IMPLEMENTATION_PLAN.md

@@ -17,6 +17,9 @@ earlier ones.
   Evidence for every change is in REVIEW.md §2.
 - 2026-09-27: Phase 10 (conformer ensembles) planned; extended the same day to
   transition structures in series and parallel channels (Dale et al. 2021).
+- 2026-09-28: Phase 10 decisions confirmed (recommended options);
+  conformer ensembles, transition structures in series, parallel channels
+  and the JSON job file implemented.
 
 **Guiding decisions**
 
@@ -458,7 +461,28 @@ a template for their own comparisons.
 **Exit criteria:** ≥ 8 reactions, every experimental number traceable to a
 DOI, REPORT.md regenerated by one command.
 
-## Phase 10 — Conformer ensembles (planned 2026-09-27)
+## Phase 10 — Conformer ensembles (planned 2026-09-27; implemented 2026-09-28)
+
+**Status (2026-09-28).** Conformer ensembles are implemented
+(`kinisot/ensemble.py`, `Conformers`, `EnsembleIsotopeEffect`,
+`equivalent_positions`; the CLI takes several files per flag with
+`--weights`, `--energies` and `--weight-uncertainty`; theory.md section
+7). Tests: `tests/test_ensemble.py`, and `tests/test_shi_ensemble.py`,
+where `compute_kie` over all 18 Shi transition structures reproduces the
+prototype to 2 × 10⁻⁶ at every site. The benchmark runner averages
+equivalent positions exactly, which moved two reported values by
+1 × 10⁻⁴ (Diels–Alder H3, nitroarene TS1B). Transition structures in
+series (`Series`) and parallel channels (`channels`) followed the same day,
+with `series_kie`/`channel_kie`, `kinisot --job` (docs/job_files.md) and
+theory.md sections 7a and 7b. Tests: `tests/test_pathways.py` (the Wittig
+and DyKAT published combinations, a three-step rate matrix to 7 × 10⁻⁸,
+channels from one reactant equal to the ensemble to 10⁻¹³) and
+`tests/test_jobs.py`. The worked example (`examples/conformers/`, made by
+`scripts/make_conformer_example.py`) followed: eight GFN2-xTB reactant
+conformers and the chair and boat transition structures, where the C4-d₂
+KIE moves by 0.010 from the lowest pair to the ensemble. Still to do: the
+Wittig and DyKAT validations from our own jobs once the logs are in
+(benchmark cases with `series` and `channels`).
 
 **Goal (requested 2026-09-27):** KIEs and EQEs from several conformers of
 each reactant, transition structure and product, given with the same atom
@@ -503,9 +527,9 @@ Four consequences shape the design:
    benchmarks). They are conformers with permuted labels and equal weights.
    - The site KIE is therefore the mean ρ over the labels divided by the
      mean ρ‡.
-   - The benchmark runner today takes the geometric mean of the pairwise
-     KIEs instead. On the Diels–Alder methyl reference, the two differ by
-     3.9 × 10⁻⁵.
+   - The benchmark runner took the geometric mean of the pairwise KIEs
+     until 2026-09-28. On the Diels–Alder methyl reference, the two differ
+     by 3.9 × 10⁻⁵.
 
 Two variants follow from the same formula:
 
@@ -745,8 +769,8 @@ dependency.
       1.033 for one oxygen and 1.066 for both, now computed from ORCA jobs.
       - In the monotopic TS1B the attacked and spectator oxygens give
         1.0474 and 1.0140. The exact average for the singly labelled
-        substrate is their harmonic mean, 1.0304; the runner's geometric
-        mean gives 1.0305.
+        substrate is their harmonic mean, 1.0304, which the runner now
+        reports (its geometric mean gave 1.0305).
       - The paper's 1.0468 is the attacked oxygen alone, which shows what
         leaving out the average costs.
     - **Channels: the DyKAT arylation** (van Dijk et al. 2021;
@@ -781,22 +805,22 @@ dependency.
       migration) are a second candidate, once their structures are
       computed.
 
-**Decisions needed**
-
-1. The default weighting: qRRHO (recommended, since conformer ensembles
-   routinely contain low modes that harmonic entropies mistreat), RRHO, or
-   requiring user free energies.
-2. The CLI syntax: `nargs="+"` per flag (recommended; backward compatible)
-   or a separate `--rct-conformers` flag.
-3. Whether the benchmarks' equivalent-position averaging moves to the exact
-   form (recommended; a 4 × 10⁻⁵ shift).
-4. Symmetry numbers: degeneracies supplied by the user only, or detected
-   with `pymsym`.
-5. Scope of the first release: parallel ensembles first, with series and
-   channels in a following release (recommended, since they reuse its
-   evaluation and result objects), or all three together.
-6. Input for series and channels on the command line: a JSON job file
-   (recommended) or new flags.
+**Decisions (confirmed 2026-09-28, the recommended options)**
+
+1. The default weighting is qRRHO, since conformer ensembles routinely
+   contain low modes that harmonic entropies mistreat; RRHO and user free
+   energies are options.
+2. The CLI takes `nargs="+"` per flag (backward compatible), not a separate
+   `--rct-conformers` flag.
+3. The benchmarks' equivalent-position averaging moves to the exact form.
+   The shifts were 1 × 10⁻⁴ at most in the reported values (Diels–Alder H3
+   and nitroarene TS1B), in their own commit with a CHANGELOG entry.
+4. Symmetry numbers and degeneracies are supplied by the user; no
+   detection.
+5. Parallel ensembles first, then series and channels reusing their
+   evaluation and result objects. Both were implemented the same day at
+   the user's request, so all three are in 2.6.0.
+6. Series and channels on the command line come from a JSON job file.
 
 **Out of scope for the first release**:
 
@@ -836,7 +860,7 @@ dependency.
 | 7 | Isotope table + syntax, projection, tunnelling, reference ✅ | v2.4.0 | ¹⁸O default change documented |
 | 8 | ASE / MLIP backend ✅ | v2.5.0 | new code path, optional extra |
 | 9 | Experimental validation suite (Singleton KIEs), scaffold ✅ | v2.5.x | needs new QC calculations and transcribed experimental values |
-| 10 | Conformer-ensemble KIEs and EQEs; transition structures in series; parallel channels | v2.7.0 | none for single-conformer inputs; equivalent-position averaging shifts by 4 × 10⁻⁵ |
+| 10 | Conformer-ensemble KIEs and EQEs, transition structures in series, parallel channels ✅ | v2.6.0 | none for single-conformer inputs; equivalent-position averaging shifted two reported values by 1 × 10⁻⁴ |
 
 Phases 0–1 are done; 2–3 are a few days and unlock adoption; 4 is the
 largest single chunk (one reviewed PR per module move); 5 is now small

README.md

@@ -79,9 +79,12 @@ The equations behind every column are in [docs/theory.md](docs/theory.md).
 ## Usage
 
 ```
-kinisot --rct FILE [--rct FILE ...] (--ts FILE | --prd FILE) --iso ATOMS [--iso ATOMS ...]
-        [-t K] [-s FACTOR] [--imag-cutoff CM-1] [--tunneling MODEL] [--barrier KCAL] [--project]
-        [--reference ATOMS] [--calc SPEC] [-o FILE] [--overwrite] [-q] [--json FILE] [--csv FILE]
+kinisot --rct FILE [FILE ...] [--rct FILE ...] (--ts FILE [FILE ...] | --prd FILE [FILE ...])
+        --iso ATOMS [--iso ATOMS ...] [-t K] [-s FACTOR] [--imag-cutoff CM-1] [--tunneling MODEL]
+        [--barrier KCAL] [--project] [--reference ATOMS] [--calc SPEC]
+        [--weights SCHEME] [--energies TABLE] [--weight-uncertainty KCAL]
+        [-o FILE] [--overwrite] [-q] [--json FILE] [--csv FILE]
+kinisot --job job.json [-o FILE] [--overwrite] [-q] [--json FILE] [--csv FILE]
 ```
 
 `python -m kinisot` is equivalent to `kinisot`.
@@ -93,6 +96,7 @@ kinisot --rct FILE [--rct FILE ...] (--ts FILE | --prd FILE) --iso ATOMS [--iso
 | KIE, one reactant | `kinisot --rct gs.out --ts ts.out --iso 5` | one `--iso` when the atom numbering is the same in both files |
 | KIE, bimolecular | `kinisot --rct dienophile.out --rct diene.out --ts ts.out --iso 0 --iso 6 --iso 15` | one `--iso` per file, `--rct` files first, `0` for a file without a substituted atom |
 | EQE | `kinisot --rct conf_a.out --prd conf_b.out --iso 24,25,26 --iso 28,29,30` | the same, with the product instead of a TS |
+| Conformers | `kinisot --rct gs_1.out gs_2.out --ts ts_*.out --iso 5` | several files after one flag are conformers of one species, with the same atom numbering; one `--iso` per flag |
 
 **Options**
 
@@ -106,8 +110,12 @@ kinisot --rct FILE [--rct FILE ...] (--ts FILE | --prd FILE) --iso ATOMS [--iso
 | `--project` | project translations and rotations out of the Hessian (Eckart) instead of discarding the six lowest modes. On by default for `--calc`/ASE inputs, off for Gaussian and ORCA (`--no-project` forces it off). |
 | `--calc SPEC`, `--delta` | compute the Hessians of geometry inputs with an ASE calculator (`xtb`, `mace_mp:medium`, `module:callable`, ...); finite-difference step in Å. |
 | `--reference ATOMS` | a second isotopologue; the KIE is also reported divided by its KIE (natural-abundance NMR style). |
+| `--weights` | how conformers are weighted: `qrrho` (default; quasi-harmonic free energies of the light isotopologue), `rrho`, `user`, `lowest` or `equal` ([theory, section 7](docs/theory.md#7-conformer-ensembles)). |
+| `--energies TABLE` | free energies of conformers, one line per file: name, G (any zero; `-` to let Kinisot compute it) and optionally a degeneracy. `--energy-unit` is `kcal/mol` (default), `kJ/mol` or `hartree`. |
+| `--weight-uncertainty` | the ensemble result gives the range of the KIE when each conformer free energy moves by this much (default 0.5 kcal/mol). |
 | `-o`, `--output` | results file (default `Kinisot_output.dat`); results are appended. `--overwrite` starts afresh, `-q` keeps the terminal quiet. |
 | `--json FILE`, `--csv FILE` | also write the full result as JSON, or append one summary row to a CSV file. |
+| `--job FILE` | a JSON job file with the structures, labels and settings: transition structures in series, parallel channels, several isotopologues in one run ([docs/job_files.md](docs/job_files.md)). The other flags then give only the outputs. |
 | `--version` | print the version. |
 
 Invalid input (an atom number out of range, a reactant with an imaginary
@@ -228,8 +236,12 @@ case: potentials whose saddle point is a different reaction.
   Wigner correction, the Skodje–Truhlar correction (which also uses the
   barrier height), or none; Bell is refused below the crossover
   temperature h c |ν‡| / 2π k.
-- Rigid rotor / harmonic oscillator, no conformational averaging, no
-  solvent corrections: compute those outside Kinisot if you need them.
+- Rigid rotor / harmonic oscillator and no solvent corrections: compute
+  those outside Kinisot if you need them.
+- Conformers given together are in fast equilibrium (Curtin–Hammett) and
+  weighted by free energies of the light isotopologue, which is exact.
+  Transition structures in series are combined at steady state, and
+  parallel channels at low conversion.
 
 ## Python API
 
@@ -263,6 +275,49 @@ subclasses); non-fatal ones are `KinisotWarning`s and are listed in
 [examples/examples.ipynb](examples/examples.ipynb). The 2.x function
 `compute_isotope_effect()` still works but is deprecated.
 
+**Conformer ensembles.** A species given as a list of files, or as
+`Conformers(files, free_energies=None, degeneracy=None)`, is an ensemble
+of conformers with the same atom numbering:
+
+```python
+from kinisot import Conformers, compute_kie
+
+r = compute_kie(rct="gs.out", ts=[["ts_chair.out", "ts_boat.out"]], iso="5", temperature=393.0)
+r.kie_tunnel                          # sum_i x_i rho_i / sum_j y_j rho'_j
+r.conformers                          # per conformer: free energy, population, its own KIE
+r.kie_tunnel_lowest, r.n_effective, r.kie_tunnel_range
+compute_kie(rct="gs.out", ts=Conformers(["a.out", "b.out"], free_energies=[0.0, 1.2]), iso="5", weights="user")
+```
+
+It returns an `EnsembleIsotopeEffect` with the same `kie`, `kie_tunnel`,
+`kie_relative`, `to_dict()` and `summary_row()`, but no ZPE/EXC/TRPF
+breakdown, which exists per conformer only. One conformer per species gives
+the ordinary `IsotopeEffect`. `kinisot.equivalent_positions(results)` gives
+the exact isotope effect for positions made equivalent by fast motion (the
+hydrogens of a rotating methyl group) from one result per placement of the
+label. [docs/theory.md, section 7](docs/theory.md#7-conformer-ensembles)
+has the equations.
+
+**Transition structures in series and parallel channels.**
+
+```python
+from kinisot import Series, channels, compute_kie, series_kie
+
+# (atom numbers illustrative) steps none of which alone commits the substrate: one label per reactant, then per step
+r = compute_kie(rct=["aldehyde.out", "ylide.out"], ts=Series(["ts_4.out", "ts_6.out"], commitment=128 / 76),
+                iso=["8", "0", "20", "20"], temperature=340.15)
+r.kie_tunnel, r.shares, r.commitment, r.commitment_for(1.033)
+# parallel routes with their own files and labels; shares from a measured selectivity
+r = channels([dict(rct="3.out", ts="ts_S.out", iso=["1", "40"]), dict(rct="3.out", ts="ts_R.out", iso=["3", "38"])],
+             shares=[3.3, 1], temperature=313.15)
+r.kie_tunnel, r.selectivity_for(1.0244)
+series_kie([1.043, 1.015], free_energies=[25.9, 26.0], temperature=340.15)   # 1.0280, from KIEs you have
+```
+
+On the command line these come from a JSON job file, `kinisot --job
+job.json` ([docs/job_files.md](docs/job_files.md)), which also runs several
+isotopologues at once.
+
 Machine-readable output from the command line: `--json run.json` (the
 full result of one run) and `--csv runs.csv` (one row per run, appended).
 
@@ -272,9 +327,12 @@ full result of one run) and `--csv runs.csv` (one row per run, appended).
   and ²H KIEs, projection, reference isotopologue, tunnelling models,
   temperature scans), a Diels–Alder reaction with one or two reactant files,
   the same KIE from ORCA and from ASE inputs (with the machine-learned
-  potential workflow), and a conformational EQE, each with the commands, the
-  expected numbers and the literature background.
+  potential workflow), a conformational EQE, and a conformer ensemble (eight
+  allyl vinyl ether conformers and the chair and boat Claisen transition
+  structures), each with the commands, the expected numbers and the
+  literature background.
 - [docs/theory.md](docs/theory.md), [docs/file_formats.md](docs/file_formats.md),
+  [docs/job_files.md](docs/job_files.md),
   [docs/faq.md](docs/faq.md), [docs/comparison.md](docs/comparison.md)
   (PyQuiver, PyQuiverHS, Gaussian's `readisotopes`, GoodVibes).
 - [benchmarks/](benchmarks/README.md): computed versus experimental KIEs,