Improvements to Chemkin file reading/parsing

rmgpy/chemkin.pyx

@@ -47,7 +47,7 @@ from rmgpy.data.kinetics.library import LibraryReaction
 from rmgpy.exceptions import ChemkinError
 from rmgpy.molecule.element import get_element
 from rmgpy.molecule.util import get_element_count
-from rmgpy.quantity import Quantity
+from rmgpy.quantity import Quantity, QuantityError
 from rmgpy.reaction import Reaction
 from rmgpy.rmg.pdep import PDepNetwork, PDepReaction
 from rmgpy.species import Species
@@ -397,7 +397,7 @@ def _read_kinetics_reaction(line, species_dict, Aunits, Aunits_surf, Eunits):
             # this identifies reactions like 'H+H+M=H2+M' as opposed to 'H+H(+M)=H2(+M)' as identified above
             third_body = True
         elif reactant not in species_dict:
-            raise ChemkinError('Unexpected reactant "{0}" in reaction {1}.'.format(reactant, reaction))
+            raise ChemkinError('Unexpected reactant "{0}" in reaction line {1}.'.format(reactant, line))
         else:
             reactant_species = species_dict[reactant]
             if not reactant_species.reactive:
@@ -1405,13 +1405,15 @@ def read_reactions_block(f, species_dict, read_comments=True):
         '{0}^3/({1}*{2})'.format(volume_units, molecule_units, time_units),  # Second-order
         '{0}^6/({1}^2*{2})'.format(volume_units, molecule_units, time_units),  # Third-order
         '{0}^9/({1}^3*{2})'.format(volume_units, molecule_units, time_units),  # Fourth-order
+        f"{volume_units}^12/({molecule_units}^4*{time_units})",  # Fifth-order
     ]
 
     Aunits_surf = [
         '',  # Zeroth-order
         's^-1'.format(time_units),  # First-order
         '{0}^2/({1}*{2})'.format(area_units, molecule_units, time_units),  # Second-order
         '{0}^4/({1}^2*{2})'.format(area_units, molecule_units, time_units),  # Third-order
+        '{0}^6/({1}^3*{2})'.format(area_units, molecule_units, time_units),  # Fourth-order
     ]
     Eunits = energy_units
 
@@ -1485,6 +1487,9 @@ def read_reactions_block(f, species_dict, read_comments=True):
                 continue
             else:
                 raise
+        except QuantityError as e:
+            logging.warning(f"Skipping the reaction {kinetics!r} due to units error: {e}")
+            continue
         reaction_list.append(reaction)
 
     return reaction_list

rmgpy/data/kinetics/family.py

@@ -2709,12 +2709,16 @@ def get_labeled_reactants_and_products(self, reactants, products, relabel_atoms=
             template_reactants = [x.item for x in template.reactants]
 
         if len(reactants0) == 1:
+            if len(template_reactants) != 1:
+                raise ValueError(f"Reaction has {len(reactants0)} reactants but template has {len(template_reactants)} reactants.")
             molecule = reactants0[0]
             mappings = self._match_reactant_to_template(molecule, template_reactants[0])
             mappings = [[map0] for map0 in mappings]
             num_mappings = len(mappings)
             reactant_structures = [molecule]
         elif len(reactants0) == 2:
+            if len(template_reactants) != 2:
+                raise ValueError(f"Reaction has {len(reactants0)} reactants but template has {len(template_reactants)} reactants.")
             molecule_a = reactants0[0]
             molecule_b = reactants0[1]
             # get mappings in forward direction
@@ -2729,6 +2733,8 @@ def get_labeled_reactants_and_products(self, reactants, products, relabel_atoms=
             reactant_structures = [molecule_a, molecule_b]
             num_mappings = len(mappings_a) * len(mappings_b)
         elif len(reactants0) == 3:
+            if len(template_reactants) != 3:
+                raise ValueError(f"Reaction has {len(reactants0)} reactants but template has {len(template_reactants)} reactants.")
             molecule_a = reactants0[0]
             molecule_b = reactants0[1]
             molecule_c = reactants0[2]
@@ -2743,7 +2749,7 @@ def get_labeled_reactants_and_products(self, reactants, products, relabel_atoms=
             reactant_structures = [molecule_a, molecule_b, molecule_c]
             num_mappings = len(mappings_a) * len(mappings_b) * len(mappings_c)
         else:
-            raise IndexError('You have {0} reactants, which is unexpected!'.format(len(reactants)))
+            raise NotImplementedError('You have {0} reactants, which is unexpected!'.format(len(reactants)))
 
         for mapping in mappings:
             try:

rmgpy/data/kinetics/rules.py

@@ -449,8 +449,8 @@ def estimate_kinetics(self, template, degeneracy=1):
         kinetics_list = remove_identical_kinetics(saved_kinetics)
 
         if len(kinetics_list) == 0:
-            raise KineticsError('Unable to determine kinetics for reaction with template {0} in family '
-                                '{1}.'.format(template, self.label))
+            raise KineticsError("Unable to determine kinetics for reaction "
+                                f"with template {template!r} in family {self.label!s}.")
 
         elif len(kinetics_list) == 1:
             kinetics, t = kinetics_list[0]

rmgpy/quantity.py

@@ -840,24 +840,23 @@ def RateCoefficient(*args, **kwargs):
 # SurfaceRateCoefficient is handled as a special case since it can take various
 # units depending on the reaction order
 SURFACERATECOEFFICIENT_CONVERSION_FACTORS = {
-    (1.0 / pq.s).dimensionality: 1.0,
-    (pq.m ** 3 / pq.s).dimensionality: 1.0,
+    (1.0 / pq.s).dimensionality: 1.0,  # unimolecular
+    (pq.m ** 3 / pq.s).dimensionality: 1.0,  # single site adsorption
     (pq.m ** 6 / pq.s).dimensionality: 1.0,
     (pq.m ** 9 / pq.s).dimensionality: 1.0,
-    (pq.m ** 3 / (pq.mol * pq.s)).dimensionality: 1.0,
+    (pq.m ** 3 / (pq.mol * pq.s)).dimensionality: 1.0,  # single site adsorption
     (pq.m ** 6 / (pq.mol ** 2 * pq.s)).dimensionality: 1.0,
     (pq.m ** 9 / (pq.mol ** 3 * pq.s)).dimensionality: 1.0,
-    (pq.m ** 2 / pq.s).dimensionality: 1.0,
-    (pq.m ** 5 / pq.s).dimensionality: 1.0,
-    (pq.m ** 2 / (pq.mol * pq.s)).dimensionality: 1.0,
-    (pq.m ** 5 / (pq.mol ** 2 * pq.s)).dimensionality: 1.0,
+    (pq.m ** 2 / pq.s).dimensionality: 1.0,  # bimolecular surface (Langmuir-Hinshelwood)
+    (pq.m ** 5 / pq.s).dimensionality: 1.0,  # dissociative adsorption
+    (pq.m ** 2 / (pq.mol * pq.s)).dimensionality: 1.0,  # bimolecular surface (Langmuir-Hinshelwood)
+    (pq.m ** 5 / (pq.mol ** 2 * pq.s)).dimensionality: 1.0,  # dissociative adsorption
     (pq.m ** 4 / (pq.mol ** 2 * pq.s)).dimensionality: 1.0,
 }
 SURFACERATECOEFFICIENT_COMMON_UNITS = [
     's^-1',  # unimolecular
     'm^3/(mol*s)', 'cm^3/(mol*s)', 'm^3/(molecule*s)', 'cm^3/(molecule*s)',  # single site adsorption
-    'm^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)',
-    # bimolecular surface (Langmuir-Hinshelwood)
+    'm^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)',  # bimolecular surface (Langmuir-Hinshelwood)
     'm^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)',  # dissociative adsorption
     'm^4/(mol^2*s)', 'cm^4/(mol^2*s)', 'm^4/(molecule^2*s)', 'cm^4/(molecule^2*s)',  # Surface_Bidentate_Dissociation
 ]