diff --git a/arkane/common.py b/arkane/common.py
index 4623a821ac..4f67ff1977 100644
--- a/arkane/common.py
+++ b/arkane/common.py
@@ -390,7 +390,11 @@ def get_element_mass(input_element, isotope=None):
         number = input_element
     elif isinstance(input_element, str):
         symbol = input_element
-        number = next(key for key, value in symbol_by_number.items() if value == input_element)
+        try:
+            number = number_by_symbol[symbol]
+        except KeyError:
+            symbol = input_element.capitalize()
+            number = number_by_symbol[symbol]
 
     if symbol is None or number is None:
         raise ValueError('Could not identify element {0}'.format(input_element))
@@ -434,6 +438,7 @@ def get_element_mass(input_element, isotope=None):
                     92: 'U', 93: 'Np', 94: 'Pu', 95: 'Am', 96: 'Cm', 97: 'Bk', 98: 'Cf', 99: 'Es', 100: 'Fm', 101: 'Md',
                     102: 'No', 103: 'Lr', 104: 'Rf', 105: 'Db', 106: 'Sg', 107: 'Bh', 108: 'Hs', 109: 'Mt', 110: 'Ds',
                     111: 'Rg', 112: 'Cn', 113: 'Nh', 114: 'Fl', 115: 'Mc', 116: 'Lv', 117: 'Ts', 118: 'Og'}
+number_by_symbol = {value: key for key, value in symbol_by_number.items()}
 
 # Structure of mass_by_symbol items: list(list(isotope1, mass1, weight1), list(isotope2, mass2, weight2), ...)
 mass_by_symbol = {
diff --git a/arkane/commonTest.py b/arkane/commonTest.py
index 7bca459d7a..50bff7fe77 100644
--- a/arkane/commonTest.py
+++ b/arkane/commonTest.py
@@ -459,6 +459,7 @@ def test_get_mass(self):
         """Test that the correct mass/number/isotope is returned from get_element_mass"""
         self.assertEqual(get_element_mass(1), (1.00782503224, 1))  # test input by integer
         self.assertEqual(get_element_mass('Si'), (27.97692653465, 14))  # test string input and most common isotope
+        self.assertEqual(get_element_mass('SI'), (27.97692653465, 14))  # test string in all caps
         self.assertEqual(get_element_mass('C', 13), (13.00335483507, 6))  # test specific isotope
         self.assertEqual(get_element_mass('Bk'), (247.0703073, 97))  # test a two-element array (no isotope data)
 
diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst
index 7b9989460d..33814e3f84 100644
--- a/documentation/source/users/rmg/input.rst
+++ b/documentation/source/users/rmg/input.rst
@@ -989,6 +989,7 @@ all of RMG's reaction families. ::
         maximumSulfurAtoms=2,
         maximumHeavyAtoms=10,
         maximumSurfaceSites=2,
+		maximumSurfaceBondOrder=4,
         maximumRadicalElectrons=2,
         maximumSingletCarbenes=1,
         maximumCarbeneRadicals=0,
diff --git a/examples/rmg/commented/input.py b/examples/rmg/commented/input.py
index 4ae3a23e6b..91b01b2b11 100644
--- a/examples/rmg/commented/input.py
+++ b/examples/rmg/commented/input.py
@@ -271,6 +271,8 @@
     maximumNitrogenAtoms=0,
     maximumSiliconAtoms=0,
     maximumSulfurAtoms=0,
+    maximumSurfaceSites=2, # maximum number of surface sites (for heterogeneous catalysis)
+    maximumSurfaceBondOrder=2, # maximum bond order of each surface sites (for heterogeneous catalysis)
     # max number of non-hydrogen atoms
     # maximumHeavyAtoms=20,
     # maximum radicals on a molecule
diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py
index 14de4275bc..bdd6607f9f 100644
--- a/rmgpy/constraints.py
+++ b/rmgpy/constraints.py
@@ -108,15 +108,21 @@ def fails_species_constraints(species):
         if struct.get_num_atoms('S') > max_sulfur_atoms:
             return True
 
+    max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1)
+    if max_heavy_atoms != -1:
+        if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms:
+            return True
+
     max_surface_sites = species_constraints.get('maximumSurfaceSites', -1)
     if max_surface_sites != -1:
         if struct.get_num_atoms('X') > max_surface_sites:
             return True
 
-    max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1)
-    if max_heavy_atoms != -1:
-        if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms:
-            return True
+    max_surface_bond_order = species_constraints.get('maximumSurfaceBondOrder', -1)
+    if max_surface_bond_order != -1:
+        for site in struct.get_surface_sites():
+            if site.get_total_bond_order() > max_surface_bond_order:
+                return True
 
     max_radicals = species_constraints.get('maximumRadicalElectrons', -1)
     if max_radicals != -1:
diff --git a/rmgpy/constraintsTest.py b/rmgpy/constraintsTest.py
index 47bceeebdb..04e6770ad9 100644
--- a/rmgpy/constraintsTest.py
+++ b/rmgpy/constraintsTest.py
@@ -62,6 +62,7 @@ def setUpClass(cls):
             maximumSiliconAtoms=1,
             maximumSulfurAtoms=1,
             maximumSurfaceSites=2,
+            maximumSurfaceBondOrder=3,
             maximumHeavyAtoms=3,
             maximumRadicalElectrons=2,
             maximumSingletCarbenes=1,
@@ -211,6 +212,16 @@ def test_surface_site_constraint(self):
         self.rmg.species_constraints['maximumCarbonAtoms'] = max_carbon
         self.rmg.species_constraints['maximumHeavyAtoms'] = max_heavy_atoms
 
+    def test_surface_bond_order_constraint(self):
+        """
+        Test that we can constrain the max bond order of surface sites.
+        """
+        mol_1site = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,Q}
+2 X u0 p0 c0 {1,Q}
+""")
+        self.assertTrue(fails_species_constraints(mol_1site))
+
     def test_heavy_constraint(self):
         """
         Test that we can constrain the max number of heavy atoms.
diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index a3d4e4b630..ff81d561c7 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -3621,9 +3621,12 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1
         inds = inds.tolist()
         revinds = [inds.index(x) for x in np.arange(len(inputs))]
 
-        pool = mp.Pool(nprocs)
+        if nprocs > 1:
+            pool = mp.Pool(nprocs)
+            kinetics_list = np.array(pool.map(_make_rule, inputs[inds]))
+        else:
+            kinetics_list = np.array(list(map(_make_rule, inputs[inds])))
 
-        kinetics_list = np.array(pool.map(_make_rule, inputs[inds]))
         kinetics_list = kinetics_list[revinds]  # fix order
 
         for i, kinetics in enumerate(kinetics_list):
@@ -4670,3 +4673,61 @@ def _child_make_tree_nodes(family, child_conn, template_rxn_map, obj, T, nprocs,
                            extension_iter_max=extension_iter_max, extension_iter_item_cap=extension_iter_item_cap)
 
     child_conn.send(list(family.groups.entries.values()))
+
+def average_kinetics(kinetics_list):
+    """
+    Based on averaging log k.
+    Hence we average n, Ea, arithmetically, but we
+    average log A (geometric average)
+    """
+    logA = 0.0
+    n = 0.0
+    Ea = 0.0
+    count = 0
+    for kinetics in kinetics_list:
+        count += 1
+        logA += np.log10(kinetics.A.value_si)
+        n += kinetics.n.value_si
+        Ea += kinetics.Ea.value_si
+
+    logA /= count
+    n /= count
+    Ea /= count
+
+    ## The above could be replaced with something like:
+    # logA, n, Ea = np.mean([[np.log10(k.A.value_si),
+    #                   k.n.value_si,
+    #                   k.Ea.value_si] for k in kinetics_list], axis=1)
+
+    Aunits = kinetics_list[0].A.units
+    if Aunits in {'cm^3/(mol*s)', 'cm^3/(molecule*s)', 'm^3/(molecule*s)'}:
+        Aunits = 'm^3/(mol*s)'
+    elif Aunits in {'cm^6/(mol^2*s)', 'cm^6/(molecule^2*s)', 'm^6/(molecule^2*s)'}:
+        Aunits = 'm^6/(mol^2*s)'
+    elif Aunits in {'s^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'}:
+        # they were already in SI
+        pass
+    elif Aunits in {'m^2/(mol*s)', 'cm^2/(mol*s)', 'm^2/(molecule*s)', 'cm^2/(molecule*s)'}:
+        # surface: bimolecular (Langmuir-Hinshelwood)
+        Aunits = 'm^2/(mol*s)'
+    elif Aunits in {'m^5/(mol^2*s)', 'cm^5/(mol^2*s)', 'm^5/(molecule^2*s)', 'cm^5/(molecule^2*s)'}:
+        # surface: dissociative adsorption
+        Aunits = 'm^5/(mol^2*s)'
+    elif Aunits == '':
+        # surface: sticking coefficient
+        pass
+    else:
+        raise Exception(f'Invalid units {Aunits} for averaging kinetics.')
+
+    if type(kinetics) not in [Arrhenius,]:
+        raise Exception(f'Invalid kinetics type {type(kinetics)!r} for {self!r}.')
+
+    if False:
+        pass
+    else:
+        averaged_kinetics = Arrhenius(
+            A=(10 ** logA, Aunits),
+            n=n,
+            Ea=(Ea * 0.001, "kJ/mol"),
+        )
+    return averaged_kinetics
diff --git a/rmgpy/data/kinetics/familyTest.py b/rmgpy/data/kinetics/familyTest.py
index af39d28659..d3b361bd40 100644
--- a/rmgpy/data/kinetics/familyTest.py
+++ b/rmgpy/data/kinetics/familyTest.py
@@ -37,12 +37,15 @@
 import numpy as np
 
 from rmgpy import settings
+import rmgpy.data.kinetics.family
 from rmgpy.data.kinetics.database import KineticsDatabase
 from rmgpy.data.kinetics.family import TemplateReaction
 from rmgpy.data.rmg import RMGDatabase
 from rmgpy.data.thermo import ThermoDatabase
 from rmgpy.molecule import Molecule
 from rmgpy.species import Species
+from rmgpy.kinetics import Arrhenius
+
 
 
 ###################################################
@@ -1090,6 +1093,23 @@ def test_retaining_atom_labels_in_template_reaction(self):
         self.assertEqual([(label, str(atom)) for label, atom in
                           reaction_list_2[0].labeled_atoms['products'].items()],
                          [('*1', 'C'), ('*2', 'C.'), ('*3', 'H')])
+        
+    def test_average_kinetics(self):
+        """
+        Test that the average kinetics are calculated correctly
+        """
+        k1 = Arrhenius(A=(1e+13, 'cm^3/(mol*s)'), n=0, Ea=(0, 'kJ/mol'), T0=(1, 'K'))
+        k2 = Arrhenius(A=(4e+13, 'cm^3/(mol*s)'), n=1, Ea=(10, 'kJ/mol'), T0=(1, 'K'))
+        kav = rmgpy.data.kinetics.family.average_kinetics([k1, k2])
+        self.assertAlmostEqual(kav.A.value_si, 2.0e7, 2) # m3/mol/s
+        self.assertAlmostEqual(kav.n.value_si, 0.5, 6)
+        self.assertAlmostEqual(kav.Ea.value_si, 5.0e3, 2)
+        self.assertAlmostEqual(kav.T0.value_si, 1.0, 6)
+        self.assertEqual(kav.A.units, 'm^3/(mol*s)')
+        self.assertAlmostEqual(np.log(kav.get_rate_coefficient(300)),
+                               np.average([np.log(k1.get_rate_coefficient(300)), 
+                                           np.log(k2.get_rate_coefficient(300))]), 6)
+
 
 
 ################################################################################
diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py
index e961e503ef..99a5e26979 100644
--- a/rmgpy/data/solvation.py
+++ b/rmgpy/data/solvation.py
@@ -1198,6 +1198,9 @@ def get_solute_data_from_groups(self, species):
         by gas-phase thermo estimate.
         """
         molecule = species.molecule[0]
+        if molecule.contains_surface_site():
+            molecule = molecule.get_desorbed_molecules()[0]
+            molecule.saturate_unfilled_valence()
         molecule.clear_labeled_atoms()
         molecule.update_atomtypes()
         solute_data = self.estimate_solute_via_group_additivity(molecule)
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index 9138b887af..271513f503 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -1374,6 +1374,7 @@ def generated_species_constraints(**kwargs):
         'maximumSiliconAtoms',
         'maximumSulfurAtoms',
         'maximumSurfaceSites',
+        'maximumSurfaceBondOrder',
         'maximumHeavyAtoms',
         'maximumRadicalElectrons',
         'maximumSingletCarbenes',
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 12c6ecdf80..b4134f60a0 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -1608,7 +1608,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False):
                                       kinetics=rxn.kinetics, duplicate=rxn.duplicate,
                                       reversible=rxn.reversible
                                       )
-            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.newReactionlist
+            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.new_reaction_list
             if getattr(r.kinetics, 'coverage_dependence', None):
                 self.process_coverage_dependence(r.kinetics)
             if not isNew:
@@ -1712,8 +1712,8 @@ def add_reaction_library_to_edge(self, reaction_library):
                                       kinetics=rxn.kinetics, duplicate=rxn.duplicate,
                                       reversible=rxn.reversible
                                       )
-            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.newReactionlist
-            if getattr(r.kinetics, 'coverage_dependence', None):
+            r, isNew = self.make_new_reaction(rxn)  # updates self.new_species_list and self.new_reaction_list
+            if r is not None and getattr(rxn.kinetics, 'coverage_dependence', None):
                 self.process_coverage_dependence(r.kinetics)
             if not isNew:
                 logging.info("This library reaction was not new: {0}".format(rxn))
@@ -1760,9 +1760,6 @@ def add_reaction_library_to_edge(self, reaction_library):
             self.add_species_to_edge(spec)
 
         for rxn in self.new_reaction_list:
-            # Note that we haven't actually evaluated any fluxes at this point
-            # Instead, we remove the comment below if the reaction is moved to
-            # the core later in the mechanism generation
             if not (self.pressure_dependence and rxn.elementary_high_p and rxn.is_unimolecular()
                     and isinstance(rxn, LibraryReaction) and isinstance(rxn.kinetics, Arrhenius) and \
                     (self.pressure_dependence.maximum_atoms is None or self.pressure_dependence.maximum_atoms >= \
diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactors.py
index 522cb33af4..d102d36c32 100644
--- a/rmgpy/rmg/reactors.py
+++ b/rmgpy/rmg/reactors.py
@@ -365,6 +365,10 @@ def __init__(self, core_phase_system, edge_phase_system, initial_conditions, ter
 
         self.terminations = [to_rms(term) for term in terminations]
 
+    def get_const_spc_indices(self,spcs=None):
+        rms_species_names = [spc.name for spc in self.core_phase_system.get_rms_species_list()]
+        return [rms_species_names.index(name) for name in self.const_spc_names]
+
     def finish_termination_criteria(self):
         """
         Convert tuples into TerminationConversion objects
diff --git a/rmgpy/yml.py b/rmgpy/yml.py
index c21e57d7d5..bd41c45d25 100644
--- a/rmgpy/yml.py
+++ b/rmgpy/yml.py
@@ -125,6 +125,7 @@ def obj_to_dict(obj, spcs, names=None, label="solvent"):
             result_dict["henrylawconstant"]["type"] = "TemperatureDependentHenryLawConstant"
             result_dict["henrylawconstant"]["Ts"] = obj.henry_law_constant_data.Ts
             result_dict["henrylawconstant"]["kHs"] = obj.henry_law_constant_data.kHs
+        result_dict["comment"] =  obj.thermo.comment
     elif isinstance(obj, NASA):
         result_dict["polys"] = [obj_to_dict(k, spcs) for k in obj.polynomials]
         result_dict["type"] = "NASA"
@@ -140,6 +141,7 @@ def obj_to_dict(obj, spcs, names=None, label="solvent"):
         result_dict["type"] = "ElementaryReaction"
         result_dict["radicalchange"] = sum([get_radicals(x) for x in obj.products]) - \
                                        sum([get_radicals(x) for x in obj.reactants])
+        result_dict["comment"] = obj.kinetics.comment
     elif isinstance(obj, Arrhenius):
         obj.change_t0(1.0)
         result_dict["type"] = "Arrhenius"