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