Merge pull request #2867 to fix constant species index mapping in SimpleReactor

Correct constant-species index mapping in SimpleReactor

Author: rwest

rmgpy/rmg/main.py

@@ -2328,7 +2328,7 @@ def calculate_cond(self, obj, Ndims, Ns=20):
         Jout /= tot  # normalize Jout
         n = self.rand_state.uniform(0, 1, 1)[0]  # draw a random number between 0 and 1
         s = 0.0
-        for indexes in np.ndenumerate(Jout):  # choose a coordinate such that grid[indexes] is choosen with probability Jout[indexes]
+        for indexes in np.ndenumerate(Jout):  # choose a coordinate such that grid[indexes] is chosen with probability Jout[indexes]
             s += Jout[indexes[0]]
             if s > n:
                 break
@@ -2402,7 +2402,7 @@ def log_conditions(rmg_memories, index):
     log newly generated reactor conditions
     """
     if rmg_memories[index].get_cond() is not None:
-        s = "conditions choosen for reactor {0} were: ".format(index)
+        s = f"conditions chosen for reactor {index} were: "
         for key, item in rmg_memories[index].get_cond().items():
             if key == "T":
                 s += "T = {0} K, ".format(item)

rmgpy/solver/simple.pyx

@@ -169,16 +169,21 @@ cdef class SimpleReactor(ReactionSystem):
                     conditions[species_dict[label]] = value
         self.sens_conditions = conditions
 
-    def get_const_spc_indices (self, core_species):
+    def get_const_spc_indices(self, core_species):
         """
         Allow to identify constant Species position in solver
         """
-        for spc in self.const_spc_names:
-            if self.const_spc_indices is None:  # initialize once the list if constant SPC declared
-                self.const_spc_indices = []
-            for spc in core_species: #Need to identify the species object corresponding to the the string written in the input file
-                if spc.label == spc:
-                    self.const_spc_indices.append(core_species.index(spc)) 
+        if self.const_spc_names is None:
+            return
+        if self.const_spc_indices is None:
+            self.const_spc_indices = []
+        else:
+            return
+        for name in self.const_spc_names:
+            for spc in core_species:
+                if spc.label == name:
+                    self.const_spc_indices.append(core_species.index(spc))
+                    break
 
     cpdef initialize_model(self, list core_species, list core_reactions, list edge_species, list edge_reactions,
                           list surface_species=None, list surface_reactions=None, list pdep_networks=None,

test/rmgpy/solver/simpleTest.py

@@ -759,3 +759,29 @@ def test_specific_collider_model(self):
         # order: Ar, N2, O2, H, CH3, CH4
         for i in range(len(simulated_mole_fracs)):
             assert round(abs(simulated_mole_fracs[i] - expected_mole_fracs[i]), 6) == 0
+
+    def test_get_const_spc_indices(self):
+        """
+        Test that const_spc_names are correctly mapped to core species indices.
+        """
+        # Dummy species with labels that should match const_spc_names
+        a = Species(smiles='C', label="CH4")
+        b = Species(smiles='[OH]', label="OH")
+        core_species = [a, b]
+
+        T = 1000.0  # K
+        P = 1.0e5   # Pa
+
+        rxn_system = SimpleReactor(
+            T,
+            P,
+            initial_mole_fractions={a: 0.5, b: 0.5},
+            n_sims=1,
+            termination=[],
+            const_spc_names=["CH4"])
+
+        # Populate const_spc_indices from labels
+        rxn_system.get_const_spc_indices(core_species)
+
+        # Only "CH4" should be marked constant
+        assert rxn_system.const_spc_indices == [core_species.index(a)]