[schneiderfelipe]

Closes #838.

Scheme stored only net stoichiometry, so species on both sides of a reaction lost their reactant coefficients. For example, 2 B -> B + C was simulated as first order in B. Scheme now stores reactants, products, transition states and equilibrium pairs, and derives everything else from them.

Fixes

• Reaction orders: get_dydt uses reactant coefficients, so 2 B -> B + C is second order in B (Wrong reaction order in an autocatalytic reaction #838).
• Molecularity: sums reactant coefficients instead of counting species, so 2 A -> B is bimolecular.
• Symmetry entropies: count molecules on each side, so 2 B -> B + C loses R ln 2.
• Repeated compounds: add up on a side, so B + B -> C is 2 B -> C (it used to be first order). Zero coefficients are dropped.
• Equilibria are paired explicitly: Scheme.equilibria stores (forward, reverse) index pairs, instead of get_k assuming consecutive half-equilibria. A -> B followed by A <=> B used to leave an unpaired half that get_k paired with the wrong reaction; now an equilibrium absorbs plain reactions between its two sides.
• Equilibrium rate constants come from equilibrium constants: get_k sets each pair from the forward half's equilibrium constant and runs Eyring only on reactions outside equilibria, so reverse halves no longer need a zero-activation special case. The resulting rate constants are unchanged.
• Symmetry entropies can't be forgotten: Scheme.reaction_freeenergies and Scheme.activation_freeenergies include them, instead of every caller subtracting temperature times symmetry entropies by hand.
• Transition states must link reactants to products: chained transition states (A -> A‡ -> A'‡ -> B), incomplete ones (A -> A‡), and transition states sharing a side or with a coefficient (C + T‡ -> D, 2 T‡) used to be silently dropped or misread; they now raise.
• Invalid schemes and inputs now raise:
  • equilibria through transition states;
  • in get_k, a reverse half of an equilibrium whose free energy is not minus the forward one, within 1e-5 J/mol;
  • negative or non-integer coefficients, and repeated compound names;
  • in the CLI and Scheme.fix, a concentration for a name that is not a compound (transition states included), a compound given both free and fixed concentrations, or a concentration without a colon; and plotting a name that is not a compound, now checked before simulating.
• Compiled models (.jk): compounds are listed in the order of Scheme.names, so models compiled by older versions stay aligned.
• Compounds are matched by name: get_k and get_kappa take values from the compounds dict by name, so a dict in any order gives the same rate constants.
• Rate constants don't depend on how free energies are given: get_k(scheme, compounds) forced every compound to the 1 M standard state, while get_freeenergies keeps gas-phase compounds at the reference pressure. Passing those free energies as delta_freeenergies gave bimolecular gas-phase rate constants 24.5 times too small (RT/p° in L/mol). get_k now converts gas-phase free energies to 1 M itself, so both paths agree, and bias now also applies with delta_freeenergies.
• Each compound keeps its own environment: get_entropies used the first compound's environment for all compounds when none was given, so mixed gas and solution models put gases at 1 M.
• Gas-phase conversions go the documented way: eyring(molecularity=) and equilibrium_constant(delta_moles=) converted from 1 M to gas-phase data, the opposite of what their docs say, off by (RT/p°)² for bimolecular reactions. get_k(volume=) was ignored; it is now the molar volume of gases.
• Unit conversions: atm-1 s-1 used RT/p instead of RT at reference pressures other than 1 atm, and unknown units passed silently for unimolecular reactions.
• Solvation models work after import overreact: solvation methods in overreact.thermo raised AttributeError unless something else had imported their module first.
• Clear errors: get_k without compounds or free energies, and solvation methods without atoms, raise TypeError instead of AssertionError. get_dydt rejects rate constants that are not one per reaction, and get_bias rejects arrays of temperatures or pressures.

Breaking changes (version bumped to 2.0.0)

Removed or made private	Use instead
parse_reactions(text)	Scheme.parse(text)
unparse_reactions(scheme)	str(scheme)
get_delta(scheme.A, x) / get_delta(scheme.B, x)	scheme.reaction_deltas(x) / scheme.activation_deltas(x)
get_reaction_entropies(transform)	scheme.reaction_entropies / scheme.activation_entropies
get_fixed_scheme(scheme, k, c)	scheme, factors = scheme.fix(c), then k * factors
get_transition_states(A, B, is_half_equilibrium)	scheme.transition_states
thermo.get_molecularity(scheme.A)	scheme.molecularity
Scheme(compounds, reactions, is_half_equilibrium, A, B)	Scheme(compounds, reactants, products, transition_states, equilibria)
scheme.B	scheme.activation_deltas
scheme.A	scheme.products - scheme.reactants
scheme.is_half_equilibrium	scheme.equilibria
is_transition_state(name)	scheme.transition_states
io.InterfaceFormatter	none (internal to the CLI)
rx.change_reference_state	rx.thermo.change_reference_state
get_kappa(scheme, compounds, method=...)	get_kappa(scheme, compounds, tunneling=...)
tunneling="none"	tunneling=None (the CLI keeps --tunneling=none)
get_bias(..., y0=[...]) in compound order	get_bias(..., y0={"A": ...}) by name
get_dydt(scheme, k, ef=...)	get_dydt(scheme, k, equilibrium_factor=...)
scheme.activation_deltas(g) - T * scheme.activation_entropies (and the reaction analogue)	scheme.activation_freeenergies(g, T) / scheme.reaction_freeenergies(g, T) (not removed, but no longer needed)

Transition states are stored apart from compounds:

• scheme.compounds no longer lists transition states, so coefficient matrices, concentration vectors (y0, get_y results) and CLI plots and tables have no transition-state entries.
• scheme.transition_states holds names (('A‡', None)) instead of indices into compounds.
• Per-name values, such as energies passed to the delta methods or bias in get_k, follow the new scheme.names: compounds, then transition states. parse_model orders model.compounds the same way, so get_freeenergies(model.compounds) stays aligned.

Scheme is now a validated frozen dataclass. Its coefficient matrices and derived data are read-only arrays, computed once and cached. Compounds are numbered in order of appearance in the resulting reactions, so Scheme.parse(str(scheme)) == scheme.

With 2-D values, the delta methods take names on the last axis, where get_delta took them on the first. get_kappa returns one row per condition instead of a flattened array.

eyring, equilibrium_constant and eckart return floats for scalar inputs, as NumPy does, instead of arrays of shape (1,) or (). The volume of equilibrium_constant and eyring is in m³/mol, as in change_reference_state.

Free energies passed to get_k follow the standard states of get_freeenergies: unsuffixed names, such as A, are gases at the reference pressure, and names such as A(w) are in solution at 1 M. Rate constants from compounds are unchanged, but those from delta_freeenergies change by a factor of RT/p°, about 24.5 L/mol, for each gas-phase molecule a step consumes, such as the second reactant of A + B -> C. To give solution free energies, name compounds with a solvent suffix.

Arrays of conditions

For example, get_k(scheme, compounds, temperature=[[300.0], [400.0]]) returns one row of rate constants per temperature. Values per compound, name or reaction lie along the last axis, and conditions broadcast against them: a scalar applies everywhere, shape (n,) gives one value per item, and shape (m, 1) adds m conditions.

• get_k takes one temperature per reaction with delta_freeenergies, or batches of shape (m, 1). It rejects temperatures that differ between compounds or between the halves of an equilibrium.
• get_internal_energies, get_enthalpies, get_entropies, get_freeenergies and get_kappa broadcast over conditions.
• The Scheme delta and free-energy methods and Scheme.fix broadcast over leading axes. The new Scheme.concentrations builds concentration arrays from a dict, by name.
• tunnel.eckart takes arrays of barriers and falls back to Wigner per element.
• get_dydt takes several states as columns, as scipy.integrate.solve_ivp does, so get_y's rate function evaluates many times in one call.
• format_numbers formats arrays with more than one dimension.
• Functions in overreact.thermo that only take scalar temperatures or pressures are documented as such.

Other changes

• Empty sides denote sources and sinks (-> A, A ->). Molecularity counts reactant molecules exactly, so sources have molecularity 0 (it used to be clamped to 1), which changes their rate constants on scales other than l mol-1 s-1.
• The CLI's activation table shows each reverse half-equilibrium's own values instead of zeros.
• Symmetry entropies and fixed concentrations are logged.
• get_dydt compiles its Jacobian once with JIT, about 100 times faster per call. get_y still doesn't pass it to the solver; re-enabling it for implicit solvers is the subject of Diffrax integration #817.
• Property-based tests with hypothesis, checked against independent oracles (invariants, round trips, detailed balance, ideal gas laws, finite differences, literature values) instead of re-deriving the code under test. Non-default tolerances are set to the measured errors.
• The coords, io, _cli and _misc modules no longer build intermediate collections or copies they don't need. Van der Waals radii, molecular volumes, the Eckart transform, Eckart quadratures and spectrum broadening use NumPy operations instead of Python loops.

Notes

• Test reliability, mypy, CI modernization, and automated releases #834: this PR touches several files that Test reliability, mypy, CI modernization, and automated releases #834 also changes; conflicts will be resolved in Test reliability, mypy, CI modernization, and automated releases #834.
• Diffrax integration #817: the Jacobian is compiled here; passing it to solvers, the 0**0 fix, trying vectorized=True in solve_ivp and an implicit default solver stay in Diffrax integration #817.
• overreact-guide: the notebooks there use parse_reactions, get_delta and get_fixed_scheme, and concentration vectors with transition-state entries, so they will need updating for 2.0.0.