[wshlavacek]

Fixes the silent half of #109. Every Arrhenius rule the load-time expansion cannot price exactly is now refused with a Tier-0 Error, which names the energy pattern or reactant context responsible, instead of running at the wrong rate. A separate PR follows with BNG-exact pricing, which lets those models run.

Why the expansion is wrong outside a narrow class

try_expand_arrhenius turns each Arrhenius rule into conventional rules with fixed rates.

• expand_binding reduces each energy pattern holding the rule's bond to at most one condition, "this site of a center molecule is bound to anything". It counts each pattern once and ORs the conditions together.
• The rule builders then build every rule from the two bond sites alone.

That is exact only when every energy pattern whose count the bond changes is one of these:

• the bond alone, A(b!1).B(a!1);
• the bond plus a lone !+ on one of its two molecules, A(b!1,c!+).B(a!1);
• the bond plus one bond to a molecule of a third type, A(b!1,c!2).B(a!1).C(a!2), where nothing but C's a can ever bind A's c.

In each case A must declare one c, not named like the reaction site. The pattern may have no state, no site required free, and no !? that BNG2 counts more than once. BNG2 divides matches by automorphisms, so k !? on m same-named components count C(m, k) times. Beyond the patterns, the rule's reactant pattern must name only the two bond sites, and a rule that breaks the bond must write + between its products.

What changed

• EnergyFunction::inexact_reason (energy_expand.cpp) checks every energy pattern against that class for the rule's reaction center.
  • A connected pattern that does not hold the bond changes its count only if it requires one of the two sites free or !+, so A(s~1), say, is left alone.
  • A pattern whose molecules are joined only by . is refused, since joining or parting any two complexes can change its count.
• collect_site_graph reads, before the rules are converted, every bond the seed species and rules' patterns can hold, and how many components of each name each molecule type declares.
• try_expand_arrhenius refuses a reactant pattern that carries anything besides the two bond sites (!? aside), and an unbinding written with .. Every refusal, old or new, now names its reason in the Error.

The nine reproducers from the issue now come back as:

case	refused with
1	energy pattern A(b!1).B(a!1,s~1) requires B.s in state 1
2	energy pattern A(b!1,c).B(a!1) requires A.c free
3	energy pattern L(r!+,r!+) constrains L.r without naming the L.r-R.l bond
4	energy pattern L(r!1,r!+).R(l!1) names a second L.r besides the reaction site
5	energy pattern L(r!1,r!2).R(l!1).S(l!2) names a second L.r besides the reaction site
6	energy pattern A(b!1,c!2).B(a!1).C(a!2) requires A.c bound to C.a, but A.c can also bind D.a
7	energy pattern A(b!1).B(a!1,c!2).C(b!2,d!3).D(c!3) names more than one bond besides the A.b-B.a bond
8	energy pattern A(b!1,c!2,d!3).B(a!1).C(a!2).D(a!3) names more than one bond besides the A.b-B.a bond
9	its reactant pattern requires A.s in state 1

A second misread in the same code: _reverse_ by name

BNG2 marks the reverse half of a reversible rule X only by naming it _reverse_X, so RM gave every Arrhenius rule of that name 1-phi. A rule the modeler labels _reverse_X is a rule as written, and BNG2 gives it phi. On bngsim's revname fixture, BNG2's network leaves 443 of 1000 complexes at t = 2; RM left 645.

A rule is now a reverse half only if the model has an Arrhenius rule X with the same rate constants and the mirror-image bond operation on the same sites, which is what BNG2 writes for a reversible rule. One case cannot be told apart from the XML: two unidirectional rules, X and _reverse_X, that mirror each other and name the same parameters for phi and Ea0. They write the same XML as one reversible rule, in either orientation, so RM reads them as one. With numbers instead of parameters, BNG2 writes per-rule _Aphi_N and _AEact0_N and the two are told apart. The documented case could only be closed by BNG2 marking reverse halves in writeXML.

Verification

• Tests.

  • arrhenius_context_test (new) loads the nine reproducers and five more shapes. Each must be refused with an Error naming the rule and the expected reason. Five models in the exact class must load with no Arrhenius diagnostic. Against main's library it fails every refusal group.

  • The five extra shapes are: a context partner only a rule makes, a context site declared twice, a !? counted twice, the other rule-context checks, and an unbinding written with ..

  • energy_expand_test covers the gate branches the fixtures do not reach.

  • The arrhenius_formula BNG2-oracle model gains seven columns, scored against a regenerated BNG2 reference (the existing columns moved in the 8th significant figure):

    • five for _reverse_: no rule X; an unrelated X; and an Arrhenius X differing from the mirror in rate constants, operation, or sites;
    • two for the lone !+.

    Main misses them by up to 102 standard errors; here the max z is 2.4.

  • ctest is 65/65 and basicmodels.py is 29/29.

• Mutation check. All 38 single-point mutations of the gate, the rule checks, the site-graph pre-pass and the reverse-half matcher are killed.

• Independent review. Two rounds by a reviewer that built its own probe models and scored them against BNG2's network.

  • Round 1 found two false accepts, !? multiplicity and _reverse_X beside an unrelated X, plus one false refusal, the lone !+. All are fixed and pinned above.
  • Round 2 found no new false accepts. It did correct the !? count from P(m, k) to C(m, k), which un-refuses L(r!1,r!?,r!?).R(l!1) on L(r,r,r) (max z 0.61 against BNG2).

• Cross-check with bngsim. bngsim's refuse_rulemonkey_energy_misreads (Re-vendor RuleMonkey 3.11.0: lift method='rm''s time() and time-TFUN refusals, run rm on absolute model time, and pick up the RuleMonkey fixes since 3.10.1 lanl/bngsim#1022) is an independent implementation of the same envelope for 3.11.0. It agrees with this gate on all 23 Arrhenius fixtures under bngsim/tests/data and on the 5 here, except for two:

  • revname, which RM now prices correctly;
  • state_change, which RM refused already.

  This gate accepts more than bngsim's in four places, each scored against BNG2's network:

  • a center type declaring the reaction site twice, max z 1.40;
  • !? counted once, !? in the rule's reactant pattern, and patterns the bond cannot change, max z 1.06, in a regime where a misapplied 0.7 energy halves the bound count;
  • the lone !+, z 1.23.

• Blast radius. 249 of the 263 XMLs under tests/ run byte-identical to main at 2 seeds, stdout and exit code. The 14 that differ are the 13 new fixtures, which main runs and this refuses, and the extended arrhenius_formula. Main's original arrhenius_formula.xml runs byte-identical under both binaries.

Once this is vendored, bngsim can drop refuse_rulemonkey_energy_misreads: RM now refuses every model that check refuses and RM misprices, with a message naming the pattern.