fix: the six minor findings of the YAML study review

apparent_activation_energy.py

@@ -209,6 +209,12 @@ def calculate_reaction_parameters(
                     else:
                         data[-1]["compound"] = key.replace(".Rate", "")
 
+        if not data:
+            raise ValueError(
+                f"no Ea could be fitted from the {'fluxes' if calculation_type == 'ri' else 'rates'}: "
+                "every one is zero or changes sign across the temperature range at the sampling "
+                "time (is the system at equilibrium by then? sample earlier)"
+            )
         return pd.DataFrame(data)
 
 

docs/USAGE.md

@@ -145,7 +145,7 @@ python microkatc.py run my_study/study.yaml       # runs every requested analysi
 python microkatc.py run my_study/study.yaml --fresh   # deletes my_study/results first
 ```
 
-`check` reports every problem at once, each with its location in the file. `run` does the same checks first.
+`check` reports every problem at once, each with its location in the file. With typed energies it also prints the barrier warnings and the overall ΔG. `run` does the same checks first.
 
 | Exit code | Meaning |
 | --- | --- |
@@ -187,7 +187,7 @@ Problems in the study file are reported before anything runs. Examples:
 | `analyses.activation_energy.sampling_time_h: 3 h is after the end of the simulation (10000 s = 2.8 h)` | Sample earlier or simulate longer. |
 | `$thermochange is not set: export thermochange=/path/to/thermochange` | Export the variable. |
 
-Warnings do not stop the run: species in no cycle (not counted in the catalyst distribution), negative barriers, and the overall ΔG from your energies (−25.9 kcal mol<sup>-1</sup> for the example, the paper's value), printed as a check.
+Warnings do not stop the run: species in no cycle (not counted in the catalyst distribution), typed energies that no step uses, negative barriers, and the overall ΔG from your energies (−25.9 kcal mol<sup>-1</sup> for the example, the paper's value), printed as a check.
 
 During a run:
 
@@ -197,6 +197,7 @@ During a run:
 | `... couldn't converge to a solution. Resimulating...` | COPASI stopped early; MicroKatc retries up to 5 times. |
 | `ConvergenceError: Solution can't be found for ...` | All retries failed: shorten the simulation or use less extreme concentrations. |
 | `C(...) = ...M couldn't reach the threshold ...` | The product never reached the conversion threshold in the simulated time; that point is left out of the conversion-time figure. Simulate longer. |
+| `no Ea could be fitted from the fluxes: ...` | Every fit was skipped, usually because the system is at equilibrium at `sampling_time_h`. Sample earlier. |
 | `Skipping Ea fit of ...: sign changes across the temperature range` | That step's net flux changes direction between temperatures, so it has no meaningful Arrhenius slope. Expected for near-equilibrium steps. |
 
 Run the test suite with `for t in tests/test_*.py; do python "$t"; done`. `tests/test_paper_barriers.py` needs `$thermochange` and is skipped without it.

microkatc.py

@@ -243,6 +243,13 @@ def main(argv=None):
     for warning in study.warnings:
         print(warning)
     if args.command == "check":
+        if study.files is None:  # typed energies: the barriers cost nothing to check
+            from thermochemistry import barrier_messages, barrier_table, gibbs_energies
+
+            gibbs = gibbs_energies(study, study.temperature_K)
+            table = barrier_table(study, gibbs)
+            for message in barrier_messages(study, gibbs, table, study.temperature_K):
+                print(message)
         print(
             f"{args.study}: OK ({len(study.steps)} steps, {len(study.species())} species)"
         )

steps.py

@@ -78,7 +78,7 @@ def parse_step(text, allow_ts=True):
     if re.search(r"\svia(\s|$)", whole):
         if not allow_ts:
             raise StepSyntaxError(f'"{whole}": via is not allowed here')
-        body, _, after = whole.partition(" via")
+        body, after = re.split(r"\svia(?:\s|$)", whole, maxsplit=1)
         names = after.split()
         if len(names) != 1:
             raise StepSyntaxError(

study.py

@@ -334,6 +334,11 @@ def load_study(path, check_energy_sources=True):
     if "temperature_K" not in conditions:
         errors.append("conditions.temperature_K: required")
     studied = conditions.get("studied_species")
+    if studied is not None and not isinstance(studied, str):
+        errors.append(
+            f"conditions.studied_species: must be one species name, got {studied!r}"
+        )
+        studied = ""
     if studied is None:
         errors.append("conditions.studied_species: required")
     elif names and studied not in names:
@@ -356,6 +361,13 @@ def load_study(path, check_energy_sources=True):
 
     # Species block
     product = species_block.get("product")
+    if product is not None and not isinstance(product, str):
+        errors.append(f"species.product: must be one species name, got {product!r}")
+        product = None
+    if product is not None and product == studied:
+        errors.append(
+            f"conditions.studied_species: {product} is the product, which always starts at 0"
+        )
     if product is not None and names and product not in names:
         errors.append(f"species.product: {product} does not appear in any step")
     overall = None
@@ -556,6 +568,11 @@ def load_study(path, check_energy_sources=True):
                 errors.append(
                     f"no energy for {name}: add it to species.energies_kcal_mol"
                 )
+        for name in energies:
+            if everything and name not in everything:
+                warnings.append(
+                    f"warning: species.energies_kcal_mol.{name} is not used by any step"
+                )
 
     # Warnings
     if cycles:

study_yaml.py

@@ -68,6 +68,10 @@ def load_yaml(path):
     try:
         with open(name, encoding="utf-8") as f:
             doc = yaml.load(f, Loader=StrictLoader)
+    except OSError as error:
+        raise StudyError(
+            [f"{name}: cannot read the file ({error.strerror})"]
+        ) from error
     except yaml.YAMLError as error:
         mark = getattr(error, "problem_mark", None) or getattr(
             error, "context_mark", None

tests/test_ea_flux_sign.py

@@ -46,6 +46,16 @@ def test_flux_sign_handling():
     assert np.isclose(fwd["Ea"], rev["Ea"]) and np.isclose(fwd["slope"], rev["slope"])
     assert np.isfinite(fwd["R2"])
 
+    nothing = df[df["name"].isin(["r03.Flux", "r04.Flux"])]
+    try:
+        ReactionParameterCalculator.calculate_reaction_parameters(
+            nothing, T, [0.1], "X", "ri", reactions=["flip", "zero"]
+        )
+    except ValueError as error:
+        assert "no Ea could be fitted" in str(error), str(error)
+    else:
+        raise AssertionError("expected ValueError when every fit is skipped")
+
     flat = np.array([1.0, 1.0, 1.0])
     assert np.isnan(PlotFunctions.calculate_R2_values(np.array(T), flat, 0.0, 1.0))
 

tests/test_microkatc.py

@@ -144,6 +144,14 @@ def broken(_study):
     assert code == 2 and "run failed: COPASI exploded" in err, err
 
 
+def test_check_shows_typed_energy_warnings():
+    study = copy_example()
+    study.write_text(study.read_text().replace("TS1: 15.0", "TS1: -5.0"))
+    code, out, _ = run("check", str(study))
+    assert code == 0 and "has a negative forward barrier (-3.0 kcal/mol)" in out, out
+    assert "Overall reaction S <=> P: ΔG = -10.0 kcal/mol at 300.0 K" in out, out
+
+
 if __name__ == "__main__":
     for name, test in list(globals().items()):
         if name.startswith("test_"):

tests/test_steps.py

@@ -70,6 +70,11 @@ def test_overall_reaction_rejects_via():
         raise AssertionError("expected StepSyntaxError")
 
 
+def test_species_name_starting_with_via():
+    step = parse_step("A + via1 <=> B via TS")
+    assert step.reactants == (("A", 1), ("via1", 1)) and step.ts == "TS", step
+
+
 if __name__ == "__main__":
     for name, test in list(globals().items()):
         if name.startswith("test_"):

tests/test_study.py

@@ -300,6 +300,33 @@ def test_reversed_and_empty_steps_are_errors():
     )
 
 
+def test_product_and_studied_species_must_be_names():
+    assert_error(
+        changed(lambda d: d["species"].update(product=["P"])),
+        "species.product: must be one species name, got ['P']",
+    )
+    assert_error(
+        changed(lambda d: d["conditions"].update(studied_species=["S"])),
+        "conditions.studied_species: must be one species name, got ['S']",
+    )
+
+
+def test_studied_species_cannot_be_the_product():
+    assert_error(
+        changed(lambda d: d["conditions"].update(studied_species="P")),
+        "conditions.studied_species: P is the product, which always starts at 0",
+    )
+
+
+def test_unused_typed_energy_is_a_warning():
+    doc = changed(lambda d: d["species"]["energies_kcal_mol"].update(TS_typo=3.0))
+    study = load_study(write(doc))
+    assert (
+        "warning: species.energies_kcal_mol.TS_typo is not used by any step"
+        in study.warnings
+    ), study.warnings
+
+
 if __name__ == "__main__":
     for name, test in list(globals().items()):
         if name.startswith("test_"):

tests/test_study_yaml.py

@@ -67,6 +67,15 @@ def test_top_level_must_be_a_mapping():
         raise AssertionError("expected StudyError")
 
 
+def test_missing_file_is_a_study_error():
+    try:
+        load_yaml(os.path.join(tempfile.mkdtemp(), "nope.yaml"))
+    except StudyError as error:
+        assert "nope.yaml: cannot read the file" in str(error), str(error)
+    else:
+        raise AssertionError("expected StudyError")
+
+
 if __name__ == "__main__":
     for name, test in list(globals().items()):
         if name.startswith("test_"):

tests/test_thermochemistry.py

@@ -205,6 +205,20 @@ def test_zero_energy_from_thermochange_stops_the_run():
     )
 
 
+def test_rerun_reports_the_same_messages():
+    doc = {
+        **TYPED,
+        "species": {
+            **TYPED["species"],
+            "energies_kcal_mol": {**TYPED["species"]["energies_kcal_mol"], "TS1": -5.0},
+        },
+    }
+    study = study_from(doc)
+    in_new_folder()
+    first = write_barrier_tables(study)
+    assert write_barrier_tables(study) == first and len(first) == 2, first
+
+
 if __name__ == "__main__":
     for name, test in list(globals().items()):
         if name.startswith("test_"):

thermochemistry.py

@@ -136,6 +136,17 @@ def overall_reaction_energy(study, gibbs):
     return products - reactants
 
 
+def barrier_messages(study, gibbs, table, temperature_K):
+    """Negative-barrier warnings, and the overall ΔG at the working temperature"""
+    messages = negative_barrier_warnings(study, table)
+    if temperature_K == study.temperature_K:
+        delta = overall_reaction_energy(study, gibbs)
+        messages.append(
+            f"Overall reaction {study.overall.text}: ΔG = {delta:.1f} kcal/mol at {temperature_K} K"
+        )
+    return messages
+
+
 def write_barrier_tables(study):
     """Writes the compound and barrier tables for every needed temperature into the current folder.
 
@@ -151,6 +162,9 @@ def write_barrier_tables(study):
             temperature
         )
         if reactions_path.is_file() and compounds_path.is_file():
+            gibbs = pd.read_csv(compounds_path, index_col=0)["Gibbs Free Energies"]
+            table = pd.read_csv(reactions_path)
+            messages += barrier_messages(study, gibbs.to_dict(), table, temperature)
             continue
         print(f"Computing Gibbs energies at {temperature} K...")
         gibbs = gibbs_energies(study, temperature)
@@ -165,10 +179,5 @@ def write_barrier_tables(study):
         partial.replace(
             reactions_path
         )  # written last and atomically: its presence means complete
-        messages += negative_barrier_warnings(study, table)
-        if temperature == study.temperature_K:
-            delta = overall_reaction_energy(study, gibbs)
-            messages.append(
-                f"Overall reaction {study.overall.text}: ΔG = {delta:.1f} kcal/mol at {temperature} K"
-            )
+        messages += barrier_messages(study, gibbs, table, temperature)
     return messages