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): """