From 4fbafa6c5878186c821d9829a5605a3093c4fd0c Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Sun, 23 Aug 2026 21:16:46 +0300 Subject: [PATCH] fix(statmech): declare reaction species in Arkane kinetics input regardless of compute_thermo render_arkane_input_template gated species(...) declarations on 'e0_only or spc.compute_thermo'. A kinetics render (e0_only=False) whose caller sets compute_thermo=False on reactants/products then emitted a reaction(...) that named species it never declared, so Arkane raised KeyError on the reactant label and produced no rate coefficient, despite all QM jobs having converged. Also declare any species named by a reaction in the render's reaction list. Thermo inputs (no reactions) are unchanged; compute_thermo still governs whether a thermo job runs, not whether a species may be declared in an input that references it. --- arc/statmech/arkane.py | 16 ++++++++++------ arc/statmech/arkane_test.py | 19 +++++++++++++++++++ 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/arc/statmech/arkane.py b/arc/statmech/arkane.py index 9f6d530ce9..6122affcbd 100644 --- a/arc/statmech/arkane.py +++ b/arc/statmech/arkane.py @@ -343,9 +343,18 @@ def render_arkane_input_template(self, skip_rotors (bool, optional): Whether to skip internal rotor consideration. Default: ``False``. e0_only (bool, optional): Whether to only run statmech (w/o thermo) to compute E0. Default: ``False``. """ + reaction_species_labels = set() + for rxn in self.reactions or list(): + reactants, products = rxn.get_reactants_and_products(return_copies=False) + rxn.reactants = [spc.label for spc in reactants] + rxn.products = [spc.label for spc in products] + reaction_species_labels.update(rxn.reactants + rxn.products) species_list = list() for spc in self.species: - if e0_only or spc.compute_thermo: + # A species named by a reaction must be declared even when its thermo is not being + # computed: Arkane's reaction() references species by label from the declared-species + # dict, so an undeclared reactant/product raises KeyError and no rate is produced. + if e0_only or spc.compute_thermo or spc.label in reaction_species_labels: smiles = spc.mol.copy(deep=True).to_smiles() if not spc.is_ts else '' adjlist = '' if smiles: @@ -387,11 +396,6 @@ def render_arkane_input_template(self, freq_scale_factor = f'\nfrequencyScaleFactor = {self.freq_scale_factor}' \ if self.freq_scale_factor is not None else '' - if self.reactions is not None: - for rxn in self.reactions: - reactants, products = rxn.get_reactants_and_products() - rxn.reactants = [spc.label for spc in reactants] - rxn.products = [spc.label for spc in products] calc_type = 'kinetics' if self.reactions else 'thermo' return Template(main_input_template).render( title=f'Arkane {calc_type} calculation', diff --git a/arc/statmech/arkane_test.py b/arc/statmech/arkane_test.py index 2afd612d23..a724f49eb4 100644 --- a/arc/statmech/arkane_test.py +++ b/arc/statmech/arkane_test.py @@ -471,6 +471,25 @@ def test_lone_pair_species_uses_adjacency_list(self): self.assertNotIn("SMILES('[CH2]')", content) # must NOT use the lossy SMILES self.assertIn("structure=SMILES('O')", content) # normal species unchanged + def test_kinetics_input_declares_reaction_species_without_compute_thermo(self): + """A reactant/product carrying ``compute_thermo=False`` must still be declared as a + ``species(...)`` line in a rendered kinetics input. Otherwise the ``reaction(...)`` block + references a species Arkane never declared, and Arkane raises ``KeyError`` on the label + (``arkane/input.py`` builds ``reactants`` from the declared-species dict). A species that no + reaction names must still be excluded by ``compute_thermo=False``.""" + rxn = ARCReaction(r_species=[ARCSpecies(label='R', smiles='C[NH]', compute_thermo=False)], + p_species=[ARCSpecies(label='P', smiles='[CH2]N', compute_thermo=False)]) + rxn.ts_species = ARCSpecies(label='TS', is_ts=True) + spc_x = ARCSpecies(label='X', smiles='O', compute_thermo=False) + adapter = ArkaneAdapter(output_directory=self.tmpdir, calcs_directory=self.tmpdir, + output_dict=dict(), sp_level=Level('gfn2'), + species=rxn.r_species + rxn.p_species + [rxn.ts_species, spc_x], + reactions=[rxn]) + content = adapter.render_arkane_input_template(statmech_dir=self.tmpdir) + self.assertIn("species('R',", content) + self.assertIn("species('P',", content) + self.assertNotIn("species('X',", content) + @classmethod def tearDownClass(cls): """