diff --git a/apparent_activation_energy.py b/apparent_activation_energy.py index acf6b20..6984f77 100644 --- a/apparent_activation_energy.py +++ b/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) diff --git a/docs/USAGE.md b/docs/USAGE.md index bb4e4b3..6b03030 100644 --- a/docs/USAGE.md +++ b/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-1 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-1 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. diff --git a/microkatc.py b/microkatc.py index 38a5bb7..c5f1eab 100644 --- a/microkatc.py +++ b/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)" ) diff --git a/steps.py b/steps.py index 119fbe8..ba07ca1 100644 --- a/steps.py +++ b/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( diff --git a/study.py b/study.py index 17c55de..db96f97 100644 --- a/study.py +++ b/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: diff --git a/study_yaml.py b/study_yaml.py index 3960c85..81d1200 100644 --- a/study_yaml.py +++ b/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 diff --git a/tests/test_ea_flux_sign.py b/tests/test_ea_flux_sign.py index 91a075f..af321f4 100644 --- a/tests/test_ea_flux_sign.py +++ b/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)) diff --git a/tests/test_microkatc.py b/tests/test_microkatc.py index 7f3ca66..e83af39 100644 --- a/tests/test_microkatc.py +++ b/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_"): diff --git a/tests/test_steps.py b/tests/test_steps.py index 77c8f6c..ce8b59f 100644 --- a/tests/test_steps.py +++ b/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_"): diff --git a/tests/test_study.py b/tests/test_study.py index 235c9f0..c984d2b 100644 --- a/tests/test_study.py +++ b/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_"): diff --git a/tests/test_study_yaml.py b/tests/test_study_yaml.py index 2393ae9..5281dd2 100644 --- a/tests/test_study_yaml.py +++ b/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_"): diff --git a/tests/test_thermochemistry.py b/tests/test_thermochemistry.py index 0482993..daac82b 100644 --- a/tests/test_thermochemistry.py +++ b/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_"): diff --git a/thermochemistry.py b/thermochemistry.py index 1d36563..dda177b 100644 --- a/thermochemistry.py +++ b/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