Correct constant-species index mapping in SimpleReactor

rmgpy/rmg/main.py

@@ -2327,7 +2327,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
@@ -2401,7 +2401,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,19 @@ 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 = []
+        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/reactionTest.py

@@ -35,6 +35,7 @@
 
 import cantera as ct
 import numpy
+import numpy as np
 import yaml
 from copy import deepcopy
 
@@ -2916,7 +2917,9 @@ def test_arrhenius(self):
             # Check that the reaction string is the same
             assert repr(converted_obj) == repr(ct_obj)
             # Check that the rate is the same. arrhenius string is not going to be identical
-            assert converted_obj.rate.input_data == ct_obj.rate.input_data
+            assert np.isclose(converted_obj.rate.input_data['rate-constant']['A'], ct_obj.rate.input_data['rate-constant']['A'])
+            assert np.isclose(converted_obj.rate.input_data['rate-constant']['b'], ct_obj.rate.input_data['rate-constant']['b'])
+            assert np.isclose(converted_obj.rate.input_data['rate-constant']['Ea'], ct_obj.rate.input_data['rate-constant']['Ea'])
 
     def test_multi_arrhenius(self):
         """
@@ -2936,9 +2939,9 @@ def test_multi_arrhenius(self):
                 # Check that the reaction string is the same
                 assert repr(converted_rxn) == repr(ct_rxn)
                 # Check that the Arrhenius rates are identical
-                assert round(abs(converted_rxn.rate.pre_exponential_factor - ct_rxn.rate.pre_exponential_factor), 3) == 0
-                assert round(abs(converted_rxn.rate.temperature_exponent - ct_rxn.rate.temperature_exponent), 7) == 0
-                assert round(abs(converted_rxn.rate.activation_energy - ct_rxn.rate.activation_energy), 7) == 0
+                assert np.isclose(converted_rxn.rate.pre_exponential_factor, ct_rxn.rate.pre_exponential_factor)
+                assert np.isclose(converted_rxn.rate.temperature_exponent, ct_rxn.rate.temperature_exponent)
+                assert np.isclose(converted_rxn.rate.activation_energy, ct_rxn.rate.activation_energy)
 
     def test_pdep_arrhenius(self):
         """
@@ -3140,18 +3143,17 @@ def test_get_reversible_potential(self):
 
     def test_get_rate_coeff(self):
         """Test get_rate_coefficient() method"""
-
         # these should be the same
         kf_1 = self.rxn_reduction.get_rate_coefficient(298,potential=0)
         kf_2 = self.rxn_reduction.kinetics.get_rate_coefficient(298,0)
 
-        assert abs(kf_1 - 43870506959779.0) < 0.000001
-        assert abs(kf_1 - kf_2) < 0.000001
+        assert np.isclose(kf_1, 43870506959779.0)
+        assert np.isclose(kf_1, kf_2)
 
-        #kf_2 should be greater than kf_1
+        # kf_2 should be greater than kf_1
         kf_1 = self.rxn_oxidation.get_rate_coefficient(298,potential=0)
         kf_2 = self.rxn_oxidation.get_rate_coefficient(298,potential=0.1)
-        assert kf_2>kf_1
+        assert kf_2 > kf_1
 
     def test_equilibrium_constant_surface_charge_transfer_kc(self):
         """

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)]