Merge pull request #2491 from rwest/matt

A collection of Matt and David's commits from #2316,
with updates and fixes added by Richard
That pull request is too large to review, so some parts are being merged first.

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 = {

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)
 

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,

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

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:

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.

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

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)
+
 
 
 ################################################################################

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)

rmgpy/rmg/input.py

@@ -1374,6 +1374,7 @@ def generated_species_constraints(**kwargs):
         'maximumSiliconAtoms',
         'maximumSulfurAtoms',
         'maximumSurfaceSites',
+        'maximumSurfaceBondOrder',
         'maximumHeavyAtoms',
         'maximumRadicalElectrons',
         'maximumSingletCarbenes',

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 >= \

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

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"