diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index de5cdf8f65..a88d4ef07b 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       max-parallel: 5
     env: # update this if needed to match a pull request on the RMG-database
-      RMG_DATABASE_BRANCH: main
+      RMG_DATABASE_BRANCH: cathub
     defaults:
       run:
         shell: bash -l {0}
diff --git a/ipython/estimate_surface_A_factors.ipynb b/ipython/estimate_surface_A_factors.ipynb
new file mode 100644
index 0000000000..d134d5a13d
--- /dev/null
+++ b/ipython/estimate_surface_A_factors.ipynb
@@ -0,0 +1,154 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "4a175662",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rmgpy.data.rmg import RMGDatabase\n",
+    "from rmgpy.kinetics import MultiArrhenius\n",
+    "from rmgpy import settings\n",
+    "\n",
+    "from surface_A_factor import *\n",
+    "from IPython.display import display\n",
+    "import numpy as np\n",
+    "import os"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "a08414a4",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "['Surface/cathub/Pd',\n",
+       " 'Surface/cathub/Ni',\n",
+       " 'Surface/cathub/Ir',\n",
+       " 'Surface/cathub/Cu',\n",
+       " 'Surface/cathub/Co',\n",
+       " 'Surface/cathub/Fe',\n",
+       " 'Surface/cathub/Ru',\n",
+       " 'Surface/cathub/Rh',\n",
+       " 'Surface/cathub/Pt',\n",
+       " 'Surface/cathub/Ag',\n",
+       " 'Surface/cathub/Au',\n",
+       " 'Surface/cathub/W']"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "libs = []\n",
+    "p = 'Surface/cathub/'\n",
+    "for e in os.listdir(os.path.join(settings['database.directory'],'kinetics','libraries',p)):\n",
+    "    libs.append(os.path.join(p,e))\n",
+    "libs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "84ef9667",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "WARNING:root:An instance of RMGDatabase already exists. Re-initializing it.\n"
+     ]
+    }
+   ],
+   "source": [
+    "database = RMGDatabase()\n",
+    "database.load(\n",
+    "    settings['database.directory'],\n",
+    "    thermo_libraries=[\n",
+    "     'primaryThermoLibrary','DFT_QCI_thermo','surfaceThermoPt111','thermo_DFT_CCSDTF12_BAC',\n",
+    "        'CHON_G4', 'NISTThermoLibrary','BurcatNS','JetSurF2.0',\n",
+    "    ],# can add more\n",
+    "    transport_libraries=[],\n",
+    "    reaction_libraries=libs,\n",
+    "    seed_mechanisms=[],\n",
+    "    kinetics_families=\"None\",\n",
+    "    kinetics_depositories=[],\n",
+    "    statmech_libraries=[],\n",
+    "    depository=[],\n",
+    "    solvation=False,\n",
+    "    surface=True,\n",
+    "    testing=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "fa5d6520",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "processing Surface/cathub/Pd...\n",
+      "processing Surface/cathub/Ni...\n",
+      "processing Surface/cathub/Ir...\n",
+      "processing Surface/cathub/Cu...\n",
+      "processing Surface/cathub/Co...\n",
+      "processing Surface/cathub/Fe...\n",
+      "processing Surface/cathub/Ru...\n",
+      "processing Surface/cathub/Rh...\n",
+      "processing Surface/cathub/Pt...\n",
+      "processing Surface/cathub/Ag...\n",
+      "processing Surface/cathub/Au...\n",
+      "processing Surface/cathub/W...\n"
+     ]
+    }
+   ],
+   "source": [
+    "for name,library in database.kinetics.libraries.items():\n",
+    "    print(f'processing {name}...')\n",
+    "    for i,entry in library.entries.items():\n",
+    "        A,comment = estimate_surface_A(entry.item)\n",
+    "        entry.long_desc += '\\n'\n",
+    "        entry.long_desc += comment\n",
+    "        if isinstance(entry.data,MultiArrhenius):\n",
+    "            for arr in entry.data.arrhenius:\n",
+    "                arr.A = A\n",
+    "        else:\n",
+    "            entry.data.A = A\n",
+    "    path = os.path.join(settings['database.directory'],'kinetics','libraries',name,'reactions.py')\n",
+    "    library.save(path)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/ipython/surface_A_factor.py b/ipython/surface_A_factor.py
new file mode 100644
index 0000000000..9f82c15315
--- /dev/null
+++ b/ipython/surface_A_factor.py
@@ -0,0 +1,172 @@
+#!/usr/bin/env python3
+
+from rmgpy.molecule import Molecule
+from rmgpy.species import Species
+from rmgpy.constants import kB,h,R
+from rmgpy.thermo.thermoengine import submit
+from rmgpy.kinetics import get_rate_coefficient_units_from_reaction_order
+
+import logging
+import numpy as np
+
+def get_A_units(rxn):
+    
+    try:
+        surf_reacts = [spcs for spcs in rxn.reactants if spcs.contains_surface_site()]
+    except IndexError:
+        surf_prods = []
+        logging.warning(f"Species do not have an rmgpy.molecule.Molecule "  
+                        "Cannot determine phases of species. We will assume gas"
+                        )
+    n_surf = len(surf_reacts)
+    n_gas = len(rxn.reactants) - len(surf_reacts)
+    return get_rate_coefficient_units_from_reaction_order(n_gas, n_surf)
+
+def get_surface_reaction_stoich(rxn):
+    
+    n_gas = 0
+    n_ads = 0
+    n_x = 0
+    
+    for s in rxn.reactants:
+        if s.is_surface_site():
+            n_x -= 1
+        elif s.contains_surface_site():
+            n_ads -= 1
+        else:
+            n_gas -= 1
+    for s in rxn.products:
+        if s.is_surface_site():
+            n_x += 1
+        elif s.contains_surface_site():
+            n_ads += 1
+        else:
+            n_gas += 1
+            
+    return (n_gas,n_ads,n_x)
+
+def get_surface_reaction_type(rxn):
+    
+    n_gas, n_ads, n_x = get_surface_reaction_stoich(rxn)
+    
+    if (n_gas, n_ads, n_x) == (1,-1,1):
+        return "desorption"
+    elif (n_gas, n_ads, n_x) == (-1,1,-1):
+        return "adsorption"
+    elif (n_gas, n_ads, n_x) == (0,1,-1):
+        return "dissociation"
+    elif (n_gas, n_ads, n_x) == (0,-1,1):
+        return "association"
+    else:
+        return "unknown"
+
+def estimate_surface_desorption_A(adsorbate: Species):
+
+    Ar_mw = 39.87750372310267 # amu
+    dummy_species = get_dummy_species(adsorbate)
+    S = dummy_species.get_entropy(298.)
+    mw = dummy_species.molecular_weight.value
+
+    A = kB*298./h
+    A *= np.exp(0.3*S/R+3.3-(18.6+np.log((mw/Ar_mw)**(5./2.)))/3)
+    commment = f"A factor estimated from gas-phase smiles {dummy_species.smiles} from "\
+    f"{dummy_species.thermo.comment} and S298={S/4.184:.2f} cal/mol/K"
+    return A,commment
+
+def get_dummy_species(adsorbate: Species):
+
+    dummy_molecules = adsorbate.molecule[0].get_desorbed_molecules()
+    for mol in dummy_molecules:
+        mol.clear_labeled_atoms()
+        
+    gas_phase_species_from_libraries = []
+    gas_phase_species_estimates = []
+    for dummy_molecule in dummy_molecules:
+        dummy_species = Species()
+        dummy_species.molecule = [dummy_molecule]
+        dummy_species.generate_resonance_structures()
+        submit(dummy_species)
+        if dummy_species.thermo.label:
+            gas_phase_species_from_libraries.append(dummy_species)
+        else:
+            gas_phase_species_estimates.append(dummy_species)
+
+    # define the comparison function to find the lowest energy
+    def lowest_energy(species):
+        if hasattr(species.thermo, 'H298'):
+            return species.thermo.H298.value_si
+        else:
+            return species.thermo.get_enthalpy(298.0)
+
+    if gas_phase_species_from_libraries:
+        species = min(gas_phase_species_from_libraries, key=lowest_energy)
+    else:
+        species = min(gas_phase_species_estimates, key=lowest_energy)
+        
+    return species
+
+def estimate_surface_dissociation_A(adsorbate: Species):
+    
+    A, comment = estimate_surface_desorption_A(adsorbate)
+    return 0.001 * A, comment
+
+def get_adsorbate(rxn, on='reactants'):
+    
+    if on == 'reactants':
+        for s in rxn.reactants:
+            if s.contains_surface_site() and not s.is_surface_site():
+                return s
+    elif on == 'products':
+        for s in rxn.products:
+            if s.contains_surface_site() and not s.is_surface_site():
+                return s
+
+def estimate_surface_A(rxn):
+    
+    reaction_type = get_surface_reaction_type(rxn)
+    comment = "A factor estimation:\n"
+    
+    if reaction_type == 'desorption':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'reactants')
+        A, _comment = estimate_surface_desorption_A(adsorbate)
+        comment += _comment
+        
+    elif reaction_type == 'dissociation':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'reactants')
+        A, _comment = estimate_surface_dissociation_A(adsorbate)
+        comment += _comment
+    
+    elif reaction_type == 'association':
+        comment += f"A factor estimate for {reaction_type}\n"
+        adsorbate = get_adsorbate(rxn,'products')
+        A, _comment = estimate_surface_dissociation_A(adsorbate)
+        comment += _comment
+        if not adsorbate.thermo:
+            submit(adsorbate)
+        deltaS = adsorbate.thermo.get_entropy(298.)
+        assert len(rxn.reactants) == 2
+        for s in rxn.reactants:
+            if s.contains_surface_site():
+                if not s.thermo:
+                    submit(s)
+                deltaS -= s.thermo.get_entropy(298.)
+        A *= np.exp(deltaS/R)
+        
+    else:
+        A = kB*298./h
+        comment = "Could not determine reaction type "\
+        f"estimating A = kb/298/h = {A:.2e}"
+        
+    units = get_A_units(rxn)
+    surface_site_density=2.483e-05 # Pt111 in metal DB
+    if units in ('m^2/(mol*s)','m^5/(mol^2*s)'):
+        A /= surface_site_density
+        comment += '\nA/=2.483e-5 mol/m^2 (Pt111 site density)'
+    elif units == 'm^4/(mol^2*s)':
+        A /= surface_site_density
+        A /= surface_site_density
+        comment += '\nA/=(2.483e-5 mol/m^2)^2 (Pt111 site density)'
+
+    return (A,units), comment
diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py
index 2b617a6b27..6001abb5f2 100644
--- a/rmgpy/data/base.py
+++ b/rmgpy/data/base.py
@@ -126,6 +126,16 @@ def __str__(self):
     def __repr__(self):
         return '<Entry index={0:d} label="{1}">'.format(self.index, self.label)
 
+    def get_metal_label(self):
+        """
+        retrieve the metal label (str) for the entry
+        """
+        if self.facet is None:
+            metal = self.metal # could be None
+        else:
+            metal = self.metal + self.facet
+        return metal
+
     def get_all_descendants(self):
         """
         retrieve all the descendants of entry
@@ -370,6 +380,12 @@ def save(self, path, reindex=True):
         f.write('#!/usr/bin/env python\n')
         f.write('# encoding: utf-8\n\n')
         f.write('name = "{0}"\n'.format(self.name))
+        if self.metal:
+            f.write('metal = "{0}"\n'.format(self.metal))
+        if self.facet:
+            f.write('facet = "{0}"\n'.format(self.facet))
+        if self.site:
+            f.write('site = "{0}"\n'.format(self.site))
         f.write('shortDesc = "{0}"\n'.format(self.short_desc))
         f.write('longDesc = """\n')
         f.write(self.long_desc.strip() + '\n')
@@ -901,8 +917,10 @@ def match_node_to_node(self, node, node_other):
         Both `node` and `node_other` must be Entry types with items containing Group or LogicNode types.
         """
         if isinstance(node.item, Group) and isinstance(node_other.item, Group):
-            return (self.match_node_to_structure(node, node_other.item, atoms=node_other.item.get_all_labeled_atoms()) and
-                    self.match_node_to_structure(node_other, node.item, atoms=node.item.get_all_labeled_atoms()))
+            return (self.match_node_to_structure(node, node_other.item, atoms=node_other.item.get_all_labeled_atoms(),
+            metal=node_other.metal, facet=node_other.facet, site=node_other.site) and
+                    self.match_node_to_structure(node_other, node.item, atoms=node.item.get_all_labeled_atoms(),
+                    metal=node.metal, facet=node.facet, site=node.site))
         elif isinstance(node.item, LogicOr) and isinstance(node_other.item, LogicOr):
             return node.item.match_logic_or(node_other.item)
         else:
@@ -918,9 +936,11 @@ def match_node_to_child(self, parent_node, child_node):
 
         if isinstance(parent_node.item, Group) and isinstance(child_node.item, Group):
             if (self.match_node_to_structure(parent_node, child_node.item,
-                                             atoms=child_node.item.get_all_labeled_atoms(), strict=True) and
+                                             atoms=child_node.item.get_all_labeled_atoms(), strict=True,
+                                             metal=child_node.metal, facet=child_node.facet, site=child_node.site) and
                     not self.match_node_to_structure(child_node, parent_node.item,
-                                                     atoms=parent_node.item.get_all_labeled_atoms(), strict=True)):
+                                                     atoms=parent_node.item.get_all_labeled_atoms(), strict=True,
+                                                     metal=parent_node.metal, facet=parent_node.facet, site=parent_node.site)):
                 return True
             return False
 
@@ -940,7 +960,7 @@ def match_node_to_child(self, parent_node, child_node):
         elif isinstance(parent_node.item, LogicOr):
             return child_node.label in parent_node.item.components
 
-    def match_node_to_structure(self, node, structure, atoms, strict=False):
+    def match_node_to_structure(self, node, structure, atoms, strict=False, metal=None, facet=None, site=None):
         """
         Return :data:`True` if the `structure` centered at `atom` matches the
         structure at `node` in the dictionary. The structure at `node` should
@@ -953,6 +973,9 @@ def match_node_to_structure(self, node, structure, atoms, strict=False):
         
         Matching to structure is more strict than to node.  All labels in structure must 
         be found in node.  However the reverse is not true, unless `strict` is set to True.
+
+        If the node has a metal attribute (metal, facet, and/or site) that does not match the
+        provdied `metal`, `facet`, and/or `site`, the structure is not a match (returns `False`)
         
         =================== ========================================================
         Attribute           Description
@@ -961,11 +984,25 @@ def match_node_to_structure(self, node, structure, atoms, strict=False):
         `structure`         A Group or a Molecule
         `atoms`             Dictionary of {label: atom} in the structure.  A possible dictionary is the one produced by structure.get_all_labeled_atoms()
         `strict`            If set to ``True``, ensures that all the node's atomLabels are matched by in the structure
+        `metal`             metal of structure (default is None)
+        `facet`             facet of structure (default is None)
+        `site`              site of structure (default is None)
         =================== ========================================================
         """
         if isinstance(node, str):
             node = self.entries[node]
         group = node.item
+        
+        if node.metal:
+            if metal != node.metal:
+                return False
+        if node.facet:
+            if facet != node.facet:
+                return False
+        if node.site:
+            if site != node.site:
+                return False
+
         if isinstance(group, LogicNode):
             return group.match_to_structure(self, structure, atoms, strict)
         else:
@@ -1023,7 +1060,7 @@ def match_node_to_structure(self, node, structure, atoms, strict=False):
 
             return result
 
-    def descend_tree(self, structure, atoms, root=None, strict=False):
+    def descend_tree(self, structure, atoms, root=None, strict=False, metal=None, facet=None, site=None):
         """
         Descend the tree in search of the functional group node that best
         matches the local structure around `atoms` in `structure`.
@@ -1035,24 +1072,28 @@ def descend_tree(self, structure, atoms, root=None, strict=False):
         Set strict to ``True`` if all labels in final matched node must match that of the
         structure.  This is used in kinetics groups to find the correct reaction template, but
         not generally used in other GAVs due to species generally not being prelabeled.
+
+        `metal` (None or str): metal of structure (default is None)
+        `facet` (None or str): facet of structure (default is None)
+        `site` (None or str): site of structure (default is None)
         """
 
         if root is None:
             for root in self.top:
-                if self.match_node_to_structure(root, structure, atoms, strict):
+                if self.match_node_to_structure(root, structure, atoms, strict, metal=metal, facet=facet, site=site):
                     break  # We've found a matching root
             else:  # didn't break - matched no top nodes
                 return None
-        elif not self.match_node_to_structure(root, structure, atoms, strict):
+        elif not self.match_node_to_structure(root, structure, atoms, strict, metal=metal, facet=facet, site=site):
             return None
 
         next_node = []
         for child in root.children:
-            if self.match_node_to_structure(child, structure, atoms, strict):
+            if self.match_node_to_structure(child, structure, atoms, strict, metal=metal, facet=facet, site=site):
                 next_node.append(child)
 
         if len(next_node) == 1:
-            return self.descend_tree(structure, atoms, next_node[0], strict)
+            return self.descend_tree(structure, atoms, next_node[0], strict, metal=metal, facet=facet, site=site)
         elif len(next_node) == 0:
             if len(root.children) > 0 and root.children[-1].label.startswith('Others-'):
                 return root.children[-1]
@@ -1062,7 +1103,7 @@ def descend_tree(self, structure, atoms, root=None, strict=False):
             # logging.warning('For {0}, a node {1} with overlapping children {2} was encountered '
             #                 'in tree with top level nodes {3}. Assuming the first match is the '
             #                 'better one.'.format(structure, root, next, self.top))
-            return self.descend_tree(structure, atoms, next_node[0], strict)
+            return self.descend_tree(structure, atoms, next_node[0], strict, metal=metal, facet=facet, site=site)
 
     def are_siblings(self, node, node_other):
         """
@@ -1144,18 +1185,22 @@ class LogicOr(LogicNode):
 
     symbol = "OR"
 
-    def match_to_structure(self, database, structure, atoms, strict=False):
+    def match_to_structure(self, database, structure, atoms, strict=False, metal=None, facet=None, site=None):
         """
         Does this node in the given database match the given structure with the labeled atoms?
         
         Setting `strict` to True makes enforces matching of atomLabels in the structure to every
         atomLabel in the node.
+
+        `metal` (None or str): metal of structure (default is None)
+        `facet` (None or str): facet of structure (default is None)
+        `site` (None or str): site of structure (default is None)
         """
         for node in self.components:
             if isinstance(node, LogicNode):
                 match = node.match_to_structure(database, structure, atoms, strict)
             else:
-                match = database.match_node_to_structure(node, structure, atoms, strict)
+                match = database.match_node_to_structure(node, structure, atoms, strict, metal=metal, facet=facet, site=site)
             if match:
                 return True != self.invert
         return False != self.invert
@@ -1195,18 +1240,22 @@ class LogicAnd(LogicNode):
 
     symbol = "AND"
 
-    def match_to_structure(self, database, structure, atoms, strict=False):
+    def match_to_structure(self, database, structure, atoms, strict=False, metal=None, facet=None, site=None):
         """
         Does this node in the given database match the given structure with the labeled atoms?
         
         Setting `strict` to True makes enforces matching of atomLabels in the structure to every
         atomLabel in the node.
+
+        `metal` (None or str): metal of structure (default is None)
+        `facet` (None or str): facet of structure (default is None)
+        `site` (None or str): site of structure (default is None)
         """
         for node in self.components:
             if isinstance(node, LogicNode):
                 match = node.match_to_structure(database, structure, atoms, strict)
             else:
-                match = database.match_node_to_structure(node, structure, atoms, strict)
+                match = database.match_node_to_structure(node, structure, atoms, strict, metal=metal, facet=facet, site=site)
             if not match:
                 return False != self.invert
         return True != self.invert
diff --git a/rmgpy/data/baseTest.py b/rmgpy/data/baseTest.py
index b1d8473d4f..b85a7e22f2 100644
--- a/rmgpy/data/baseTest.py
+++ b/rmgpy/data/baseTest.py
@@ -61,6 +61,17 @@ def test_match_node_to_structure(self):
                 """)
         )
 
+        entry1_facet = Entry(
+            item=Group().from_adjacency_list(
+                """
+                1 *3 C  u1 {2,D} {3,S}
+                2    C  u0 {1,D}
+                3 *5 Cd u0 {1,S} {4,D}
+                4    C  u0 {3,D}
+                """),
+            facet='111'
+        )
+
         entry2 = Entry(
             item=Group().from_adjacency_list(
                 """
@@ -92,6 +103,17 @@ def test_match_node_to_structure(self):
         # entry3 contains fewer labels than entry1, therefore it can be matched
         self.assertTrue(self.database.match_node_to_structure(entry1, entry3.item, atoms=entry3.item.get_all_labeled_atoms()))
 
+        # entry3 does not have a facet, so it should not match entry1_facet
+        self.assertFalse(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms()))
+
+        # entry3 should match entry_1 facet if we provide a matching facet
+        self.assertTrue(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms(),
+        facet='111'))
+
+        # entry3 should not match entry1_facet if the facets do not match
+        self.assertFalse(self.database.match_node_to_structure(entry1_facet, entry3.item, atoms=entry3.item.get_all_labeled_atoms(),
+        facet='211'))
+
     def test_match_node_to_node(self):
         """
         Test that nodes can match other nodes.
@@ -103,14 +125,25 @@ def test_match_node_to_node(self):
                 """)
         )
 
+        entry1_facet = Entry(
+            item=Group().from_adjacency_list(
+                """
+                1 *1 R!H u1
+                """),
+            facet = '111'
+        )
+
         entry2 = Entry(
             item=Group().from_adjacency_list(
                 """
                 1 *1 Cb u1
                 """)
         )
+
         self.assertTrue(self.database.match_node_to_node(entry1, entry1))
+        self.assertTrue(self.database.match_node_to_node(entry1_facet, entry1_facet))
         self.assertFalse(self.database.match_node_to_node(entry1, entry2))
+        self.assertFalse(self.database.match_node_to_node(entry1, entry1_facet))
 
 
 class TestForbiddenStructures(unittest.TestCase):
diff --git a/rmgpy/data/kinetics/database.py b/rmgpy/data/kinetics/database.py
index d37e0de432..aebb83568a 100644
--- a/rmgpy/data/kinetics/database.py
+++ b/rmgpy/data/kinetics/database.py
@@ -44,7 +44,8 @@
 from rmgpy.kinetics import Arrhenius, ArrheniusEP, ThirdBody, Lindemann, Troe, \
                            PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \
                            Chebyshev, KineticsData, StickingCoefficient, \
-                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, ArrheniusBM
+                           StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP, \
+                           ArrheniusBM, SurfaceArrheniusBM
 from rmgpy.molecule import Molecule, Group
 from rmgpy.reaction import Reaction, same_species_lists
 from rmgpy.species import Species
@@ -79,9 +80,13 @@ def __init__(self):
             'SurfaceArrhenius': SurfaceArrhenius,
             'SurfaceArrheniusBEP': SurfaceArrheniusBEP,
             'R': constants.R,
-            'ArrheniusBM': ArrheniusBM
+            'ArrheniusBM': ArrheniusBM,
+            'SurfaceArrheniusBM': SurfaceArrheniusBM
         }
         self.global_context = {}
+        self.metal = None
+        self.facet = None
+        self.site = None
 
     def __reduce__(self):
         """
@@ -558,7 +563,9 @@ def react_molecules(self, molecules, products=None, only_families=None, prod_res
             if only_families is None or label in only_families:
                 try:
                     reaction_list.extend(family.generate_reactions(molecules, products=products,
-                                                                   prod_resonance=prod_resonance))
+                                                                   prod_resonance=prod_resonance,
+                                                                   metal=self.metal, facet=self.facet,
+                                                                   site=self.site))
                 except:
                     logging.error("Problem family: {}".format(label))
                     logging.error("Problem reactants: {}".format(molecules))
@@ -813,6 +820,6 @@ def reconstruct_kinetics_from_source(self, reaction, source, fix_barrier_height=
             else:
 
                 h298 = rxn_copy.get_enthalpy_of_reaction(298)
-                if isinstance(kinetics, (ArrheniusEP, ArrheniusBM)):
+                if isinstance(kinetics, (ArrheniusEP, ArrheniusBM, SurfaceArrheniusBM)):
                     kinetics = kinetics.to_arrhenius(h298)
                 return kinetics
diff --git a/rmgpy/data/kinetics/depository.py b/rmgpy/data/kinetics/depository.py
index ad24c7db36..32263ce310 100644
--- a/rmgpy/data/kinetics/depository.py
+++ b/rmgpy/data/kinetics/depository.py
@@ -221,7 +221,8 @@ def load_entry(self,
         """
 
         reaction = Reaction(reactants=[], products=[], specific_collider=specificCollider,
-                            degeneracy=degeneracy, duplicate=duplicate, reversible=reversible)
+                            degeneracy=degeneracy, duplicate=duplicate, reversible=reversible, metal=metal,
+                            facet=facet, site=site)
 
         entry = Entry(
             index=index,
diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 29794d7875..d884a9a665 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -55,7 +55,7 @@
 from rmgpy.exceptions import ActionError, DatabaseError, InvalidActionError, KekulizationError, KineticsError, \
                              ForbiddenStructureException, UndeterminableKineticsError
 from rmgpy.kinetics import Arrhenius, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficient, \
-                           StickingCoefficientBEP, ArrheniusBM
+                           StickingCoefficientBEP, ArrheniusBM, SurfaceArrheniusBM
 from rmgpy.kinetics.uncertainties import RateUncertainty, rank_accuracy_map
 from rmgpy.molecule import Bond, GroupBond, Group, Molecule
 from rmgpy.molecule.atomtype import ATOMTYPES
@@ -727,16 +727,19 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
         self.boundary_atoms = local_context.get('boundaryAtoms', None)
         self.tree_distances = local_context.get('treeDistances', None)
         self.reverse_map = local_context.get('reverseMap', None)
+        self.vdw_bonds = local_context.get('vdwBonds', None)
 
         self.reactant_num = local_context.get('reactantNum', None)
         self.product_num = local_context.get('productNum', None)
 
         self.auto_generated = local_context.get('autoGenerated', False)
 
-        if self.reactant_num:
-            self.groups.reactant_num = self.reactant_num
-        else:
-            self.groups.reactant_num = len(self.forward_template.reactants)
+        if not self.reactant_num:
+            self.reactant_num = len(self.forward_template.reactants)
+        if not self.product_num:
+            self.product_num = len(self.forward_template.products)
+        
+        self.groups.reactant_num = self.reactant_num
 
         # Generate the reverse template if necessary
         self.forward_template.reactants = [self.groups.entries[label] for label in self.forward_template.reactants]
@@ -858,8 +861,37 @@ def save_entry(self, f, entry):
         """
         return save_entry(f, entry)
 
+    def split_template(self):
+        """
+        If the family has one template and is not unimolecular,
+        split the reaction template into multiple reactants
+        """
+        if self.reactant_num > 1 and len(self.forward_template.reactants) == 1:
+            template = self.forward_template.reactants[0]
+            if isinstance(template.item, Group):
+                if self.vdw_bonds:
+                    for label0,label1 in self.vdw_bonds.items():
+                        atom0 = template.item.get_labeled_atoms(label0)[0]
+                        atom1 = template.item.get_labeled_atoms(label1)[0]
+                        bond = GroupBond(atom0, atom1, order=[0])
+                        template.item.add_bond(bond)
+
+                groups = template.item.split()
+                template_reactants = []
+                if len(groups) > 1:
+                    for i,group in enumerate(groups):
+                        group.remove_van_der_waals_bonds()
+                        template_reactant = deepcopy(template)
+                        template_reactant.item = group
+                        template_reactant.label += str(i)
+                        template_reactants.append(template_reactant)
+                self.forward_template.reactants = template_reactants
+                if self.reverse_template:
+                    self.reverse_template.products = template_reactants
+                assert len(self.forward_template.reactants) == self.reactant_num
+
     def save_training_reactions(self, reactions, reference=None, reference_type='', short_desc='', long_desc='',
-                                rank=3):
+                                metal=None, facet=None, site=None, rank=3):
         """
         This function takes a list of reactions appends it to the training reactions file.  It ignores the existence of
         duplicate reactions.  
@@ -878,6 +910,12 @@ def save_training_reactions(self, reactions, reference=None, reference_type='',
             short_desc = [short_desc] * len(reactions)
         if not isinstance(long_desc, list):
             long_desc = [long_desc] * len(reactions)
+        if not isinstance(metal, list):
+            metal = [metal] * len(reactions)
+        if not isinstance(facet, list):
+            facet = [facet] * len(reactions)
+        if not isinstance(site, list):
+            site = [site] * len(reactions)
         if not isinstance(rank, list):
             rank = [rank] * len(reactions)
 
@@ -948,6 +986,9 @@ def save_training_reactions(self, reactions, reference=None, reference_type='',
                 short_desc=str(short_desc[i]),
                 long_desc=str(long_desc[i]),
                 rank=rank[i],
+                metal=metal[i],
+                facet=facet[i],
+                site=site[i]
             )
 
             # Add this entry to the loaded depository so it is immediately usable
@@ -1014,6 +1055,9 @@ def save_groups(self, path):
         if self.reverse_map is not None:
             f.write('reverseMap = {0}\n\n'.format(self.reverse_map))
 
+        if self.vdw_bonds is not None:
+            f.write('vdwBonds = {0}\n\n'.format(self.vdw_bonds))
+
         if self.reactant_num is not None:
             f.write('reactantNum = {0}\n\n'.format(self.reactant_num))
         if self.product_num is not None:
@@ -1272,10 +1316,7 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
                 if quantum_mechanics:
                     quantum_mechanics.run_jobs(item.reactants + item.products, procnum=procnum)
 
-            if entry.facet is None:
-                metal = entry.metal # could be None
-            else:
-                metal = entry.metal + entry.facet
+            metal = entry.get_metal_label()
 
             for reactant in item.reactants:
                 # Clear atom labels to avoid effects on thermo generation, ok because this is a deepcopy
@@ -1812,7 +1853,8 @@ def _match_reactant_to_template(self, reactant, template_reactant):
         else:
             raise NotImplementedError("Not expecting template of type {}".format(type(struct)))
 
-    def generate_reactions(self, reactants, products=None, prod_resonance=True, delete_labels=True, relabel_atoms=True):
+    def generate_reactions(self, reactants, products=None, prod_resonance=True, delete_labels=True, relabel_atoms=True,
+                           metal= None, facet=None, site=None):
         """
         Generate all reactions between the provided list of one, two, or three
         `reactants`, which should be either single :class:`Molecule` objects
@@ -1847,6 +1889,9 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele
                                      prod_resonance=prod_resonance,
                                      delete_labels=delete_labels,
                                      relabel_atoms=relabel_atoms,
+                                     metal=metal,
+                                     facet=facet,
+                                     site=site
                                      ))
 
         if not self.own_reverse and self.reversible:
@@ -1858,6 +1903,9 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele
                                          prod_resonance=prod_resonance,
                                          delete_labels=delete_labels,
                                          relabel_atoms=relabel_atoms,
+                                         metal=metal,
+                                         facet=facet,
+                                         site=site
                                          ))
         return reaction_list
 
@@ -1869,6 +1917,9 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True):
         Returns `True` if successful and `False` if the reverse reaction is forbidden.
         Will raise a `KineticsError` if unsuccessful for other reasons.
         """
+
+        #parse out metal attrs of rxn
+        metal, facet, site = rxn.metal, rxn.facet, rxn.site
         if self.own_reverse and all([spc.has_reactive_molecule() for spc in rxn.products]):
             # Check if the reactants are the same
             same_reactants = 0
@@ -1888,7 +1939,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True):
 
             reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products],
                                                      products=rxn.reactants, forward=True,
-                                                     react_non_reactive=react_non_reactive)
+                                                     react_non_reactive=react_non_reactive,
+                                                     metal=metal, facet=facet, site=site)
             reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self)
             if len(reactions) == 0:
                 logging.error("Expecting one matching reverse reaction, not zero in reaction family {0} for "
@@ -1911,7 +1963,8 @@ def add_reverse_attribute(self, rxn, react_non_reactive=True):
                 try:
                     reaction_list = self._generate_reactions([spc.molecule for spc in rxn.products],
                                                              products=rxn.reactants, forward=True,
-                                                             react_non_reactive=react_non_reactive)
+                                                             react_non_reactive=react_non_reactive,
+                                                             metal=metal, facet=facet, site=site)
                     reactions = find_degenerate_reactions(reaction_list, same_reactants, kinetics_family=self)
                 finally:
                     self.forbidden = temp_object
@@ -1960,6 +2013,8 @@ def calculate_degeneracy(self, reaction):
         identical reactants (since you will be reducing them later in the algorithm), add
         `ignoreSameReactants= True` to this method.
         """
+        # parse out metal attrs of reaction
+        metal, facet, site = reaction.metal, reaction.facet, reaction.site
         # Check if the reactants are the same
         # If they refer to the same memory address, then make a deep copy so
         # they can be manipulated independently
@@ -2004,7 +2059,7 @@ def calculate_degeneracy(self, reaction):
         reactions = []
         for combo in molecule_combos:
             reactions.extend(self._generate_reactions(combo, products=reaction.products, forward=True,
-                                                      react_non_reactive=True))
+                                                      react_non_reactive=True, metal=metal, facet=facet, site=site))
 
         # remove degenerate reactions
         reactions = find_degenerate_reactions(reactions, same_reactants, template=reaction.template,
@@ -2022,7 +2077,8 @@ def calculate_degeneracy(self, reaction):
         return reactions[0].degeneracy
 
     def _generate_reactions(self, reactants, products=None, forward=True, prod_resonance=True,
-                            react_non_reactive=False, delete_labels=True, relabel_atoms=True):
+                            react_non_reactive=False, delete_labels=True, relabel_atoms=True,
+                            metal=None, facet=None, site=None):
         """
         Generate a list of all the possible reactions of this family between
         the list of `reactants`. The number of reactants provided must match
@@ -2060,6 +2116,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
         # This also makes a copy of the reactants list so we don't modify the
         # original
         reactants = [reactant if isinstance(reactant, list) else [reactant] for reactant in reactants]
+        self.split_template()
 
         if forward:
             template = self.forward_template
@@ -2073,21 +2130,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
         if self.auto_generated and reactant_num != len(reactants):
             return []
 
-        if len(reactants) > len(template.reactants):
-            # If the template contains a surface site, we do not want to split it because it will break vdw bonds
-            if isinstance(template.reactants[0].item, Group):
-                if template.reactants[0].item.contains_surface_site():
-                    return []
-            # if the family has one template and is bimolecular split template into multiple reactants
-            try:
-                grps = template.reactants[0].item.split()
-                template_reactants = []
-                for grp in grps:
-                    template_reactants.append(grp)
-            except AttributeError:
-                template_reactants = [x.item for x in template.reactants]
-        else:
-            template_reactants = [x.item for x in template.reactants]
+        template_reactants = [x.item for x in template.reactants]
 
         # Unimolecular reactants: A --> products
         if len(reactants) == 1 and len(template_reactants) == 1:
@@ -2409,6 +2452,12 @@ def generate_products_and_reactions(order):
         # Also store the reaction template (useful so we can easily get the kinetics later)
         for reaction in rxn_list:
 
+            if any([metal,facet,site]):
+                if reaction.is_surface_reaction():
+                    reaction.metal = metal
+                    reaction.facet = facet
+                    reaction.site = site
+
             # Restore the labeled atoms long enough to generate some metadata
             for reactant in reaction.reactants:
                 reactant.clear_labeled_atoms()
@@ -2653,19 +2702,29 @@ def get_kinetics_from_depository(self, depository, reaction, template, degenerac
                     '[{0}]'.format(';'.join([g.label for g in template])))
                 kinetics.comment += "\nfamily: {}".format(self.label)
                 if entry.metal:
-                    kinetics.comment += "\nmetal: {}".format(self.metal)
+                    kinetics.comment += "\nmetal: {}".format(entry.metal)
                 if entry.facet:
-                    kinetics.comment += "\nfacet: {}".format(self.facet)
-                if entry.facet:
-                    kinetics.comment += "\nsite: {}".format(self.site)
+                    kinetics.comment += "\nfacet: {}".format(entry.facet)
+                if entry.site:
+                    kinetics.comment += "\nsite: {}".format(entry.site)
         return kinetics_list
 
-    def _select_best_kinetics(self, kinetics_list):
+    def _select_best_kinetics(self, kinetics_list, metal=None, facet=None):
         """
         For a given set of kinetics `kinetics_list`, return the kinetics deemed
         to be the "best". This is determined to be the one with the lowest
         non-zero rank that occurs first (has the lowest index).
         """
+        if metal:
+            # select only the reactions that match metals
+            _kinetics_list = [k for k in kinetics_list if k[1].metal == metal]
+            if _kinetics_list:
+                kinetics_list = _kinetics_list
+        if facet:
+            # select only the reactions that match facets
+            _kinetics_list = [k for k in kinetics_list if k[1].facet == facet]
+            if _kinetics_list:
+                kinetics_list = _kinetics_list
         if any([x[1].rank == 0 for x in kinetics_list]) and not all([x[1].rank == 0 for x in kinetics_list]):
             kinetics_list = [x for x in kinetics_list if x[1].rank != 0]
         kinetics_list.sort(key=lambda x: (x[1].rank, x[1].index))
@@ -2695,11 +2754,14 @@ def get_kinetics(self, reaction, template_labels, degeneracy=1, estimator='', re
 
         template = self.retrieve_template(template_labels)
 
+        metal = reaction.metal
+        facet = reaction.facet
+
         # Check the various depositories for kinetics
         for depository in depositories:
             kinetics_list0 = self.get_kinetics_from_depository(depository, reaction, template, degeneracy)
             if len(kinetics_list0) > 0 and not return_all_kinetics:
-                kinetics, entry, is_forward = self._select_best_kinetics(kinetics_list0)
+                kinetics, entry, is_forward = self._select_best_kinetics(kinetics_list0, metal=metal, facet=facet)
                 return kinetics, depository, entry, is_forward
             else:
                 for kinetics, entry, is_forward in kinetics_list0:
@@ -2814,23 +2876,16 @@ def get_labeled_reactants_and_products(self, reactants, products):
         new entities in memory so input molecules `reactants` and `products` won't be affected.
         If RMG cannot find appropriate labels, (None, None) will be returned.
         """
+
+        self.split_template() # if the family has one template and is bimolecular split template into multiple reactants
         template = self.forward_template
         reactants0 = [reactant.copy(deep=True) for reactant in reactants]
+        template_reactants = [x.item for x in template.reactants]
 
         if self.auto_generated and self.reactant_num != len(reactants):
             return None, None
-
-        if len(reactants) > len(template.reactants):
-            # if the family has one template and is bimolecular split template into multiple reactants
-            try:
-                grps = template.reactants[0].item.split()
-                template_reactants = []
-                for grp in grps:
-                    template_reactants.append(grp)
-            except AttributeError:
-                template_reactants = [x.item for x in template.reactants]
-        else:
-            template_reactants = [x.item for x in template.reactants]
+        
+        template_reactants = [x.item for x in template.reactants]
 
         if len(reactants0) == 1:
             molecule = reactants0[0]
@@ -2969,14 +3024,14 @@ def get_training_depository(self):
         else:
             raise DatabaseError('Could not find training depository in family {0}.'.format(self.label))
 
-    def add_entry(self, parent, grp, name):
+    def add_entry(self, parent, grp, name, facet=None):
         """
         Adds a group entry with parent parent
         group structure grp
         and group name name
         """
         ind = len(self.groups.entries) - 1
-        entry = Entry(index=ind, label=name, item=grp, parent=parent)
+        entry = Entry(index=ind, label=name, item=grp, parent=parent, facet=facet)
         self.groups.entries[name] = entry
         self.rules.entries[name] = []
         if entry.parent:
@@ -3217,9 +3272,12 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np.
 
         if exts == [] and not gave_up_split:  # should only occur when all reactions at this node are identical
             rs = template_rxn_map[parent.label]
+            split_violation_flag = False
             for q, rxn in enumerate(rs):
+                if split_violation_flag:
+                    break
                 for j in range(q):
-                    if not same_species_lists(rxn.reactants, rs[j].reactants, generate_initial_map=True, save_order=self.save_order):
+                    if (not split_violation_flag) and (not same_species_lists(rxn.reactants, rs[j].reactants, generate_initial_map=True, save_order=self.save_order)):
                         for p, atm in enumerate(parent.item.atoms):
                             if atm.reg_dim_atm[0] != atm.reg_dim_atm[1]:
                                 logging.error('atom violation')
@@ -3262,11 +3320,32 @@ def extend_node(self, parent, template_rxn_map, obj=None, T=1000.0, iter_max=np.
                             for prod in rxn.products:
                                 logging.error(prod.label)
                                 logging.error(prod.to_adjacency_list())
-                        logging.error("Clearing Regularization Dimensions and Reattempting")  # this usually happens when node expansion breaks some symmetry
-                        parent.item.clear_reg_dims()  # this almost always solves the problem
-                        return True
+                        if 'split_violations' in parent.item.props:
+                            parent.item.props['split_violations'] += 1
+                        else:
+                            parent.item.props['split_violations'] = 1
+                        split_violation_flag = True
+                        if parent.item.props['split_violations'] >= 5:
+                            # give up
+                            # break the loop, and add facets to extend the tree
+                            logging.error(f"Exceeded max split violations for {parent.label}")
+                            break
+                        else:
+                            logging.error("Clearing Regularization Dimensions and Reattempting")  # this usually happens when node expansion breaks some symmetry
+                            parent.item.clear_reg_dims()  # this almost always solves the problem
+                            return True
+            # we either have all matching reactions or we broke the loop and exceeded max split violations for parent node
+            # lets try to extend the tree by splitting based on facet
+            if not parent.facet:
+                facets = list(set(r.facet for r in rs if r.facet))
+                if len(facets) > 1: # we have more than 1 facet to split
+                    for facet in facets:
+                        # if (not parent.facet) or (parent.facet == facet):
+                        #     # the parent does not have a facet or the facets match, so lets create a child node with the facet
+                            group = parent.item
+                            extname = parent.label + '_facet{}'.format(facet)
+                            self.add_entry(parent, group, extname, facet=facet)
             return False
-        
         if gave_up_split:
             return False
         
@@ -3670,9 +3749,8 @@ 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)
-
-        kinetics_list = np.array(pool.map(_make_rule, inputs[inds]))
+        with mp.Pool(nprocs) as pool:
+            kinetics_list = np.array(pool.map(_make_rule, inputs[inds]))
         kinetics_list = kinetics_list[revinds]  # fix order
 
         for i, kinetics in enumerate(kinetics_list):
@@ -3708,6 +3786,11 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1
         if template_rxn_map is None:
             template_rxn_map = self.get_reaction_matches(remove_degeneracy=True, get_reverse=True, fix_labels=True)
 
+        if template_rxn_map['Root'][0].is_surface_reaction():
+            surface = True
+        else:
+            surface = False
+
         rxns = np.array(template_rxn_map['Root'])
 
         if test_rxn_inds is None:
@@ -3759,7 +3842,10 @@ def cross_validate(self, folds=5, template_rxn_map=None, test_rxn_inds=None, T=1
                     L = list(set(template_rxn_map[entry.label]) - set(rxns_test))
 
                     if L != []:
-                        kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions)
+                        if surface:
+                            kinetics = SurfaceArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions)
+                        else:
+                            kinetics = ArrheniusBM().fit_to_reactions(L, recipe=self.forward_recipe.actions)
                         kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(T))
                         k = kinetics.get_rate_coefficient(T)
                         errors[rxn] = np.log(k / krxn)
@@ -3940,7 +4026,7 @@ def simple_regularization(self, node, template_rxn_map, test=True):
                                 if not self.rxns_match_node(node, rxns):
                                     atm1.radical_electrons = oldvals
 
-                if (not skip and atm1.reg_dim_r[1] != [] and
+                if (not grp.contains_surface_site() and not skip and atm1.reg_dim_r[1] != [] and
                         ('inRing' not in atm1.props.keys() or atm1.reg_dim_r[1][0] != atm1.props['inRing'])):
                     if 'inRing' not in atm1.props.keys():
                         if (all(['inRing' in child.item.atoms[i].props.keys() for child in node.children]) and
@@ -4162,19 +4248,25 @@ def get_reactant_thermo(reactant,metal):
         for i, entry in enumerate(entries):
             if estimate_thermo:
                 # parse out the metal to scale to
-                if entry.facet is None:
-                    metal = entry.metal # could be None
-                else:
-                    metal = entry.metal + entry.facet
+                metal = entry.get_metal_label()
                 for j, react in enumerate(entry.item.reactants):
+                    if metal:
+                        # we need to make a new species for different metals beacause thermo is different
+                        rxns[i].reactants[j] = react.copy(deep=True)
+                        rxns[i].reactants[j].thermo = None
                     if rxns[i].reactants[j].thermo is None:
-                        rxns[i].reactants[j].thermo = get_reactant_thermo(react,metal)
+                        rxns[i].reactants[j].thermo = get_reactant_thermo(rxns[i].reactants[j],metal)
 
                 for j, react in enumerate(entry.item.products):
+                    if metal:
+                        # we need to make a new species for different metals beacause thermo is different
+                        rxns[i].products[j] = react.copy(deep=True)
+                        rxns[i].products[j].thermo = None
                     if rxns[i].products[j].thermo is None:
-                        rxns[i].products[j].thermo = get_reactant_thermo(react,metal)
+                        rxns[i].products[j].thermo = get_reactant_thermo(rxns[i].products[j],metal)
             rxns[i].kinetics = entry.data
             rxns[i].rank = entry.rank
+            rxns[i].index = entry.index
 
             if remove_degeneracy:  # adjust for degeneracy
                 rxns[i].kinetics.A.value_si /= rxns[i].degeneracy
@@ -4245,13 +4337,18 @@ def get_reactant_thermo(reactant,metal):
                     reacts = [Species(molecule=[get_label_fixed_mol(x.molecule[0], root_labels)], thermo=x.thermo)
                               for x in rxns[i].reactants]
                     rrev = Reaction(reactants=products, products=reacts,
-                                    kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank)
+                                    kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank,
+                                    metal=rxns[i].metal, facet=rxns[i].facet, site=rxns[i].site)
                     rrev.is_forward = False
 
                     if estimate_thermo:
-                        for rev_react in rrev.reactants:
-                            if rev_react.thermo is None:
-                                rev_react.thermo = get_reactant_thermo(rev_react,metal)
+                        for x,react in enumerate(rrev.reactants):
+                            if metal:
+                            # we need to make a new species for different metals beacause thermo is different
+                                rrev.reactants[x] = react.copy(deep=True)
+                                rrev.reactants[x].thermo = None
+                            if rrev.reactants[x].thermo is None:
+                                rrev.reactants[x].thermo = get_reactant_thermo(rrev.reactants[x],metal)
 
                     rev_rxns.append(rrev)
 
@@ -4283,15 +4380,21 @@ def get_reactant_thermo(reactant,metal):
                     products = [Species(molecule=[p]) for p in products]
 
                 rrev = Reaction(reactants=products, products=rxns[i].reactants,
-                                kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank)
+                                kinetics=rxns[i].generate_reverse_rate_coefficient(), rank=rxns[i].rank,
+                                metal = rxns[i].metal, facet = rxns[i].facet, site = rxns[i].site)
 
                 rrev.is_forward = False
 
                 if estimate_thermo:
-                    for rev_react in rrev.reactants:
-                        if rev_react.thermo is None:
-                            rev_react.thermo = get_reactant_thermo(rev_react,metal)
+                    for x,react in enumerate(rrev.reactants):
+                        if metal:
+                        # we need to make a new species for different metals beacause thermo is different
+                            rrev.reactants[x] = react.copy(deep=True)
+                            rrev.reactants[x].thermo = None
+                        if rrev.reactants[x].thermo is None:
+                            rrev.reactants[x].thermo = get_reactant_thermo(rrev.reactants[x],metal)
                 rxns[i] = rrev
+                rxns[i].index = entry.index
 
         if self.own_reverse and get_reverse:
             return rxns + rev_rxns
@@ -4325,9 +4428,9 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac
                 else:
                     mol = mol.merge(r.molecule[0])
             try:
-                flag = not self.is_entry_match(mol, root, resonance=True)
+                flag = not self.is_entry_match(mol, root, resonance=True, metal=rxn.metal, facet=rxn.facet, site=rxn.site)
             except:
-                flag = not self.is_entry_match(mol, root, resonance=False)
+                flag = not self.is_entry_match(mol, root, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site)
 
             if flag:
                 logging.error(root.item.to_adjacency_list())
@@ -4344,7 +4447,7 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac
 
             while entry.children != []:
                 for child in entry.children:
-                    if self.is_entry_match(mol, child, resonance=False):
+                    if self.is_entry_match(mol, child, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site):
                         entry = child
                         rxn_lists[child.label].append(rxn)
                         break
@@ -4362,10 +4465,21 @@ def get_reaction_matches(self, rxns=None, thermo_database=None, remove_degenerac
 
         return rxn_lists
 
-    def is_entry_match(self, mol, entry, resonance=True):
+    def is_entry_match(self, mol, entry, resonance=True, metal=None, facet=None, site=None):
         """
         determines if the labeled molecule object of reactants matches the entry entry
         """
+
+        if entry.metal:
+            if metal != entry.metal:
+                return False
+        if entry.facet:
+            if facet != entry.facet:
+                return False
+        if entry.site:
+            if site != entry.site:
+                return False
+
         if isinstance(entry.item, Group):
             if resonance:
                 structs = mol.generate_resonance_structures()
@@ -4385,7 +4499,7 @@ def rxns_match_node(self, node, rxns):
                 else:
                     mol = mol.merge(r.molecule[0])
 
-            if not self.is_entry_match(mol, node, resonance=False):
+            if not self.is_entry_match(mol, node, resonance=False, metal=rxn.metal, facet=rxn.facet, site=rxn.site):
                 return False
 
         return True
@@ -4663,14 +4777,31 @@ def _make_rule(rr):
     """
     recipe, rxns, Tref, fmax, label, ranks = rr
     n = len(rxns)
+    if n == 0:
+        return None
     for i, rxn in enumerate(rxns):
         rxn.rank = ranks[i]
+    if rxn.is_surface_reaction():
+        surface = True
+    else:
+        surface = False
     rxns = np.array(rxns)
     data_mean = np.mean(np.log([r.kinetics.get_rate_coefficient(Tref) for r in rxns]))
-    if n > 0:
+    if surface:
+        kin = SurfaceArrheniusBM().fit_to_reactions(rxns, recipe=recipe)
+    else:
         kin = ArrheniusBM().fit_to_reactions(rxns, recipe=recipe)
-        if n == 1:
-            kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label)
+    if n == 1:
+        kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label)
+    else:
+        if surface:
+            dlnks = np.array([
+                np.log(
+                    SurfaceArrheniusBM().fit_to_reactions(rxns[list(set(range(len(rxns))) - {i})], recipe=recipe)
+                    .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref))
+                    .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)
+                ) for i, rxn in enumerate(rxns)
+            ])  # 1) fit to set of reactions without the current reaction (k)  2) compute log(kfit/kactual) at Tref
         else:
             dlnks = np.array([
                 np.log(
@@ -4679,17 +4810,15 @@ def _make_rule(rr):
                     .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)
                 ) for i, rxn in enumerate(rxns)
             ])  # 1) fit to set of reactions without the current reaction (k)  2) compute log(kfit/kactual) at Tref
-            varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2
-            # weighted average calculations
-            ws = 1.0 / varis
-            V1 = ws.sum()
-            V2 = (ws ** 2).sum()
-            mu = np.dot(ws, dlnks) / V1
-            s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1))
-            kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label)
-        return kin
-    else:
-        return None
+        varis = (np.array([rank_accuracy_map[rxn.rank].value_si for rxn in rxns]) / (2.0 * 8.314 * Tref)) ** 2
+        # weighted average calculations
+        ws = 1.0 / varis
+        V1 = ws.sum()
+        V2 = (ws ** 2).sum()
+        mu = np.dot(ws, dlnks) / V1
+        s = np.sqrt(np.dot(ws, (dlnks - mu) ** 2) / (V1 - V2 / V1))
+        kin.uncertainty = RateUncertainty(mu=mu, var=s ** 2, N=n, Tref=Tref, data_mean=data_mean, correlation=label)
+    return kin
 
 
 def _spawn_tree_process(family, template_rxn_map, obj, T, nprocs, depth, min_splitable_entry_num, min_rxns_to_spawn, extension_iter_max, extension_iter_item_cap):
diff --git a/rmgpy/data/kinetics/groups.py b/rmgpy/data/kinetics/groups.py
index 55c47bd263..405bc4dbd0 100644
--- a/rmgpy/data/kinetics/groups.py
+++ b/rmgpy/data/kinetics/groups.py
@@ -76,7 +76,7 @@ def __repr__(self):
         return '<KineticsGroups "{0}">'.format(self.label)
 
     def load_entry(self, index, label, group, kinetics, reference=None, referenceType='', shortDesc='', longDesc='',
-                   nodalDistance=None):
+                   nodalDistance=None, metal=None, facet=None, site=None):
         """
         Method for parsing entries in database files.
         Note that these argument names are retained for backward compatibility.
@@ -103,7 +103,10 @@ def load_entry(self, index, label, group, kinetics, reference=None, referenceTyp
             reference_type=referenceType,
             short_desc=shortDesc,
             long_desc=longDesc.strip(),
-            nodal_distance=nodalDistance
+            nodal_distance=nodalDistance,
+            metal = metal,
+            facet = facet,
+            site = site
         )
 
     def get_reaction_template(self, reaction):
@@ -115,6 +118,7 @@ def get_reaction_template(self, reaction):
 
         # Get forward reaction template and remove any duplicates
         forward_template = self.top[:]
+        metal, facet, site = reaction.metal, reaction.facet, reaction.site
 
         temporary = []
         symmetric_tree = False
@@ -148,7 +152,7 @@ def get_reaction_template(self, reaction):
 
             atoms = r.get_all_labeled_atoms()
 
-            matched_node = self.descend_tree(r, atoms, root=entry, strict=True)
+            matched_node = self.descend_tree(r, atoms, root=entry, strict=True, metal=metal, facet=facet, site=site)
 
             if matched_node is not None:
                 template.append(matched_node)
@@ -175,7 +179,7 @@ def get_reaction_template(self, reaction):
                     # Match structures
                     atoms = reactant.get_all_labeled_atoms()
                     # Descend the tree, making sure to match atomlabels exactly using strict = True
-                    matched_node = self.descend_tree(reactant, atoms, root=entry, strict=True)
+                    matched_node = self.descend_tree(reactant, atoms, root=entry, strict=True, metal=metal, facet=facet, site=site)
                     if matched_node is not None:
                         template.append(matched_node)
                     # else:
diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py
index 220b813994..fc59c772ac 100644
--- a/rmgpy/data/kinetics/library.py
+++ b/rmgpy/data/kinetics/library.py
@@ -295,11 +295,16 @@ def check_for_duplicates(self, mark_duplicates=False):
         """
         for entry0 in self.entries.values():
             reaction0 = entry0.item
+            metal0 = entry0.get_metal_label()
             if not reaction0.duplicate:
                 # This reaction is not marked as a duplicate reaction
                 # This means that if we find any duplicate reactions, it is an error
                 for entry in self.entries.values():
                     reaction = entry.item
+                    metal = entry.get_metal_label()
+                    if metal0 != metal:
+                        # the metals don't match, so they are not duplicate reactions
+                        continue
                     if reaction0 is not reaction and reaction0.is_isomorphic(reaction):
                         # We found a duplicate reaction that wasn't marked!
                         # RMG requires all duplicate reactions to be marked, unlike CHEMKIN
@@ -545,7 +550,7 @@ def load_entry(self,
         # Make a blank reaction
         rxn = Reaction(reactants=[], products=[], degeneracy=degeneracy, duplicate=duplicate, reversible=reversible,
                        allow_pdep_route=allow_pdep_route, elementary_high_p=elementary_high_p,
-                       allow_max_rate_violation=allow_max_rate_violation)
+                       allow_max_rate_violation=allow_max_rate_violation,metal=metal, facet=facet, site=site)
         # if not rxn.is_balanced():
         #    raise DatabaseError('Reaction {0} in kinetics library {1} was not balanced! Please reformulate.'.format(rxn, self.label))
         # label = str(rxn)
diff --git a/rmgpy/data/surface.py b/rmgpy/data/surface.py
index a4559adaed..cf52d1662d 100644
--- a/rmgpy/data/surface.py
+++ b/rmgpy/data/surface.py
@@ -247,9 +247,9 @@ def find_binding_energies(self, metal):
                 metal_binding_energies = self.libraries['surface'].get_binding_energies(metal)
             except DatabaseError:
                 # no exact match was found, so continue on as if no facet was given
-                logging.warning("Requested metal %r not found in database.", metal)
+                logging.info("Requested metal %r not found in database.", metal)
                 metal = metal[:facet.span()[0]]
-                logging.warning("Searching for generic %r.", metal)
+                logging.info("Searching for generic %r.", metal)
                 facet = None
 
         if facet is None:
@@ -263,7 +263,7 @@ def find_binding_energies(self, metal):
             else:  # multiple matches
                 # average the binding energies together? just pick the first one?
                 # just picking the first one for now...
-                logging.warning(f"Found multiple binding energies for {metal!r}. Using {metal_entry_matches[0]!r}.")
+                logging.info(f"Found multiple binding energies for {metal!r}. Using {metal_entry_matches[0]!r}.")
                 metal_binding_energies = self.libraries['surface'].get_binding_energies(metal_entry_matches[0])
 
         return metal_binding_energies
diff --git a/rmgpy/kinetics/__init__.py b/rmgpy/kinetics/__init__.py
index a9dc32e49b..10bc31ef48 100644
--- a/rmgpy/kinetics/__init__.py
+++ b/rmgpy/kinetics/__init__.py
@@ -35,4 +35,4 @@
 from rmgpy.kinetics.kineticsdata import KineticsData, PDepKineticsData
 from rmgpy.kinetics.tunneling import Wigner, Eckart
 from rmgpy.kinetics.surface import SurfaceArrhenius, SurfaceArrheniusBEP, \
-                     StickingCoefficient, StickingCoefficientBEP
+                    StickingCoefficient, StickingCoefficientBEP, SurfaceArrheniusBM
diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx
index 244318a61d..7190f13154 100644
--- a/rmgpy/kinetics/arrhenius.pyx
+++ b/rmgpy/kinetics/arrhenius.pyx
@@ -34,8 +34,9 @@ from scipy.optimize import curve_fit, fsolve
 
 cimport rmgpy.constants as constants
 import rmgpy.quantity as quantity
-from rmgpy.exceptions import KineticsError
+from rmgpy.exceptions import KineticsError, InvalidActionError
 from rmgpy.kinetics.uncertainties import rank_accuracy_map
+from rmgpy.kinetics import get_rate_coefficient_units_from_reaction_order
 from rmgpy.molecule.molecule import Bond
 
 # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value
@@ -669,9 +670,11 @@ cdef class ArrheniusBM(KineticsModel):
             self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts)
 
         # fill in parameters
-        A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)']
-        order = len(rxns[0].reactants)
-        self.A = (A, A_units[order])
+        surf_reacts = [spcs for spcs in rxns[0].reactants if spcs.contains_surface_site()]
+        n_surf = len(surf_reacts)
+        n_gas = len(rxns[0].reactants) - len(surf_reacts)
+        A_units = get_rate_coefficient_units_from_reaction_order(n_gas, n_surf)
+        self.A = (A, A_units)
 
         self.n = n
         self.w0 = (w0, 'J/mol')
@@ -1106,6 +1109,8 @@ def get_w0(actions, rxn):
 
     recipe = actions
 
+    metal = rxn.get_metal_label()
+
     wb = 0.0
     wf = 0.0
     for act in recipe:
@@ -1120,14 +1125,24 @@ def get_w0(actions, rxn):
                 atom2 = a_dict[act[3]]
 
         if act[0] == 'BREAK_BOND':
-            bd = mol.get_bond(atom1, atom2)
-            wb += bd.get_bde()
+            bd = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
+            wb += bd.get_bde(metal=metal)
         elif act[0] == 'FORM_BOND':
-            bd = Bond(atom1, atom2, act[2])
-            wf += bd.get_bde()
+            bd = Bond(a_dict[act[1]], a_dict[act[3]], act[2])
+            wf += bd.get_bde(metal=metal)
         elif act[0] == 'CHANGE_BOND':
-            bd1 = mol.get_bond(atom1, atom2)
-
+            if not mol.has_bond(a_dict[act[1]], a_dict[act[3]]): # we dont have a bond
+                is_vdW_bond = False
+                for atom in (a_dict[act[1]], a_dict[act[3]]):
+                    if atom.is_surface_site():
+                        is_vdW_bond = True
+                        break
+                if not is_vdW_bond: # no surface site, so no vdW bond
+                    raise InvalidActionError('Attempted to change a nonexistent bond.')
+                else: # we found a surface site, so we will make vdw bond
+                    bd1 = Bond(a_dict[act[1]], a_dict[act[3]], order=0)
+            else: # we have a bond
+                bd1 = mol.get_bond(a_dict[act[1]], a_dict[act[3]])
             if act[2] + bd1.order == 0.5:
                 mol2 = None
                 for r in rxn.products:
@@ -1151,8 +1166,9 @@ def get_w0(actions, rxn):
             if bd2.order == 0:
                 bd2_bde = 0.0
             else:
-                bd2_bde = bd2.get_bde()
-            bde_diff = bd2_bde - bd1.get_bde()
+                bd2_bde = bd2.get_bde(metal=metal)
+
+            bde_diff = bd2_bde - bd1.get_bde(metal=metal)
             if bde_diff > 0:
                 wf += abs(bde_diff)
             else:
diff --git a/rmgpy/kinetics/surface.pxd b/rmgpy/kinetics/surface.pxd
index 8e37ac54b7..bbd87ca5b9 100644
--- a/rmgpy/kinetics/surface.pxd
+++ b/rmgpy/kinetics/surface.pxd
@@ -28,7 +28,7 @@
 cimport numpy as np
 
 from rmgpy.kinetics.model cimport KineticsModel
-from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP
+from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP, ArrheniusBM
 from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity
 
 ################################################################################
@@ -72,4 +72,7 @@ cdef class SurfaceArrhenius(Arrhenius):
 cdef class SurfaceArrheniusBEP(ArrheniusEP):
     cdef public dict _coverage_dependence
     pass
+################################################################################
+cdef class SurfaceArrheniusBM(ArrheniusBM):
+    pass
 
diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx
index 2c108e9204..3bf48cd4f0 100644
--- a/rmgpy/kinetics/surface.pyx
+++ b/rmgpy/kinetics/surface.pyx
@@ -625,3 +625,76 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP):
             coverage_dependence=self.coverage_dependence,
             comment=self.comment,
         )
+
+################################################################################
+
+cdef class SurfaceArrheniusBM(ArrheniusBM):
+    """
+    A kinetics model based on the (modified) Arrhenius equation, using the
+    Blowers-Masel equation to determine the activation energy.
+    Based on Blowers and Masel's 2000 paper Engineering Approximations for Activation
+    Energies in Hydrogen Transfer Reactions.
+    The attributes are:
+
+    =============== =============================================================
+    Attribute       Description
+    =============== =============================================================
+    `A`             The preexponential factor
+    `n`             The temperature exponent
+    `w0`            The average of the bond dissociation energies of the bond formed and the bond broken
+    `E0`            The activation energy for a thermoneutral reaction
+    `Tmin`          The minimum temperature at which the model is valid, or zero if unknown or undefined
+    `Tmax`          The maximum temperature at which the model is valid, or zero if unknown or undefined
+    `Pmin`          The minimum pressure at which the model is valid, or zero if unknown or undefined
+    `Pmax`          The maximum pressure at which the model is valid, or zero if unknown or undefined
+    `comment`       Information about the model (e.g. its source)
+    =============== =============================================================
+
+    """
+
+    property A:
+        """The preexponential factor. 
+    
+        This is the only thing different from a ArrheniusBM class"""
+        def __get__(self):
+            return self._A
+        def __set__(self, value):
+            self._A = quantity.SurfaceRateCoefficient(value)
+
+    def __repr__(self):
+        """
+        Return a string representation that can be used to reconstruct the
+        ArrheniusBM object.
+        """
+        string = 'SurfaceArrheniusBM(A={0!r}, n={1!r}, w0={2!r}, E0={3!r}'.format(self.A, self.n, self.w0, self.E0)
+        if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin)
+        if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax)
+        if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin)
+        if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax)
+        if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty)
+        if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
+        string += ')'
+        return string
+
+    def __reduce__(self):
+        """
+        A helper function used when pickling an ArrheniusEP object.
+        """
+        return (SurfaceArrheniusBM, (self.A, self.n, self.w0, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax,
+                              self.uncertainty, self.comment))
+
+    cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn):
+        """
+        Return an :class:`Arrhenius` instance of the kinetics model using the
+        given enthalpy of reaction `dHrxn` to determine the activation energy.
+        """
+        return SurfaceArrhenius(
+            A=self.A,
+            n=self.n,
+            Ea=(self.get_activation_energy(dHrxn) * 0.001, "kJ/mol"),
+            T0=(1, "K"),
+            Tmin=self.Tmin,
+            Tmax=self.Tmax,
+            uncertainty=self.uncertainty,
+            comment=self.comment,
+        )
\ No newline at end of file
diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index cb720ef84f..ede13fa00e 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -120,6 +120,8 @@ cdef class Bond(Edge):
 
     cpdef bint is_van_der_waals(self) except -2
 
+    cpdef bint is_surface_bond(self) except -2
+
     cpdef bint is_single(self) except -2
 
     cpdef bint is_double(self) except -2
@@ -170,6 +172,8 @@ cdef class Molecule(Graph):
 
     cpdef remove_van_der_waals_bonds(self)
 
+    cpdef add_van_der_waals_bond(self)
+
     cpdef sort_atoms(self)
     
     cpdef str get_formula(self)
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index d0a057bc2f..1e98a15e01 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -45,13 +45,14 @@
 import cython
 import numpy as np
 
+import rmgpy.quantity as quantity
 import rmgpy.constants as constants
 import rmgpy.molecule.converter as converter
 import rmgpy.molecule.element as elements
 import rmgpy.molecule.group as gr
 import rmgpy.molecule.resonance as resonance
 import rmgpy.molecule.translator as translator
-from rmgpy.exceptions import DependencyError
+from rmgpy.exceptions import DependencyError, DatabaseError
 from rmgpy.molecule.adjlist import Saturator
 from rmgpy.molecule.atomtype import AtomType, ATOMTYPES, get_atomtype, AtomTypeError
 from rmgpy.molecule.element import bdes
@@ -638,16 +639,74 @@ def sorting_key(self):
                 self.atom1.number if self.atom1 is not None else 0,
                 self.atom2.number if self.atom2 is not None else 0)
 
-    def get_bde(self):
+    def get_bde(self, metal='Pt111'):
         """
         estimate the bond dissociation energy in J/mol of the bond based on the order of the bond
         and the atoms involved in the bond
+
+        If the bond involves a surface site, the atomic binding energies from the `metal` (str) 
+        will be use to estimate the bond dissociation energy.
         """
         try:
             return bdes[(self.atom1.element.symbol, self.atom2.element.symbol, self.order)]
         except KeyError:
-            raise KeyError('Bond Dissociation energy not known for combination: '
-                           '({0},{1},{2})'.format(self.atom1.element.symbol, self.atom2.element.symbol, self.order))
+            if not self.is_surface_bond():
+                raise KeyError('Bond Dissociation energy not known for combination: '
+                               '({0},{1},{2})'.format(self.atom1.element.symbol, self.atom2.element.symbol, self.order))
+            else:
+                if self.atom1.is_surface_site():
+                    adatom = self.atom2.element.symbol
+                else:
+                    adatom = self.atom1.element.symbol
+                if self.is_van_der_waals():
+                    if adatom == 'O':
+                        binding_energy = quantity.Energy(-0.189,'eV/molecule') # binding energy for H2O_ads in surfaceThermoPt111
+                    elif adatom == 'N':
+                        binding_energy = quantity.Energy(-0.673,'eV/molecule') # binding energy for NH3_ads in surfaceThermoPt111
+                    elif adatom == 'C':
+                        binding_energy = quantity.Energy(-0.122,'eV/molecule') # binding energy for CH4_ads in surfaceThermoPt111
+                    elif adatom == 'H':
+                        binding_energy = quantity.Energy(-0.054,'eV/molecule') # binding energy for H2_ads in surfaceThermoPt111
+                    else:
+                        binding_energy = quantity.Energy(-0.2,'eV/molecule') # default vdw binding energy
+                    return -1 * binding_energy.value_si
+                try:
+                    from rmgpy.data.rmg import get_db
+                    tdb = get_db('thermo')
+                    metal_db = tdb.surface['metal']
+                except (DatabaseError, KeyError):
+                    from rmgpy.data.surface import MetalDatabase
+                    from rmgpy import settings
+                    metal_db = MetalDatabase()
+                    metal_db.load(os.path.join(settings['database.directory'], 'surface'))
+                metal_binding_energies = metal_db.find_binding_energies(metal=metal)
+                binding_energy = metal_binding_energies[adatom]
+                max_bond_order = {'C': 4., 'O': 2., 'N': 3., 'H': 1.}
+                normalized_bond_order = self.order / max_bond_order[adatom]
+                # intercepts estimated from DOI:10.1103/PhysRevLett.99.016105
+                intercepts = {
+                    'C' : {
+                        1 : 0,
+                        2 : -1.3,
+                        3 : -1.1,
+                        4 : 0,
+                    },
+                    'N' : {
+                        1 : -0.5,
+                        2 : -0.75,
+                        3 : 0
+                    },
+                    'O' : {
+                        1 : -0.5,
+                        2 : 0
+                    },
+                    'H' : {
+                        1 : 0
+                    }
+                }
+                intercept = quantity.Energy(intercepts[adatom][self.order],'eV/molecule')
+                binding_energy = quantity.Energy(binding_energy).value_si * normalized_bond_order + intercept.value_si
+                return -1 * binding_energy
 
     def equivalent(self, other):
         """
@@ -752,6 +811,16 @@ def is_van_der_waals(self):
         """
         return self.is_order(0)
 
+    def is_surface_bond(self):
+        """
+        Return ``True`` if the bond represents a surface bond or 
+        ``False`` if not.
+        """
+        if self.atom1.is_surface_site() or self.atom2.is_surface_site():
+            return True
+        else:
+            return False
+
     def is_order(self, other_order):
         """
         Return ``True`` if the bond is of order other_order or ``False`` if
@@ -1129,6 +1198,45 @@ def remove_van_der_waals_bonds(self):
             if bond.is_van_der_waals():
                 self.remove_bond(bond)
 
+    def add_van_der_waals_bond(self):
+        """
+        Adds a single van der Waals bond to the graph to connect a surface site to the physisorbed absorbate.
+        If the adsorbate has lone pairs, the vdw bond will be added to an atom with lone pairs.
+        Atoms will be selected in the folowing order:
+            1) Nitrogen with lone pairs
+            2) Oxygen with lone pairs
+            3) atoms without lone pairs (C and H)
+        """
+        cython.declare(bond=Bond, fragment=Molecule, adsorbate=Molecule,
+            surface_site=Molecule,fragments=list, has_surface_site=cython.bint, atom=Atom, lone_pair_atoms=list)
+
+        if self.contains_surface_site():
+            fragments = self.split()
+            if len(fragments) == 2:
+                has_surface_site = False
+                for fragment in fragments:
+                    if fragment.is_surface_site():
+                        has_surface_site = True
+                        surface_site = fragment
+                    else:
+                        adsorbate = fragment
+                if has_surface_site:
+                    lone_pair_atoms = [atom for atom in adsorbate.atoms if atom.lone_pairs > 0]
+                    if len(lone_pair_atoms) == 1:
+                        # only one atom with a lone pair, so pick this one
+                        atom = lone_pair_atoms[0]
+                    elif len(lone_pair_atoms) > 1:
+                        # more than one atom with lone pair
+                        # sort by element symbol and pick the first (pick nitrogen before oxygen)
+                        lone_pair_atoms.sort(key = lambda atom: atom.symbol)
+                        atom = lone_pair_atoms[0]
+                    else:
+                        # no lone pairs, so pick the first atom
+                        atom = adsorbate.atoms[0]
+                    # create and add the vdw bond
+                    bond = Bond(atom, surface_site.atoms[0],0.0)
+                    self.add_bond(bond)
+
     def sort_atoms(self):
         """
         Sort the atoms in the graph. This can make certain operations, e.g.
diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py
index 576d73c374..0372e51688 100644
--- a/rmgpy/molecule/moleculeTest.py
+++ b/rmgpy/molecule/moleculeTest.py
@@ -840,6 +840,20 @@ def test_get_bond_string(self):
         bond = Bond(atom1=Atom(element=get_element(1)), atom2=Atom(element=get_element(6)), order=1)
         self.assertEqual(bond.get_bond_string(), 'C-H')
 
+    def test_get_bde(self):
+
+        bonds = [Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=0),
+                Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=1),
+                Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=2),
+                Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=3),
+                Bond(atom1=Atom(element=get_element('X')), atom2=Atom(element=get_element(6)), order=4)]
+
+        bde = 0
+        for bond in bonds:
+            bde_bond = bond.get_bde('Pt111')
+            self.assertGreater(bde_bond, bde)
+            bde = bde_bond
+
 ################################################################################
 
 
@@ -2759,6 +2773,16 @@ def test_remove_van_der_waals_bonds(self):
         mol.remove_van_der_waals_bonds()
         self.assertEqual(len(mol.get_all_edges()), 1)
 
+    def test_add_van_der_waals_bond(self):
+        """Test we can add a van-der-Waals bond"""
+        mol = Molecule(smiles='*.NOC')
+        mol.add_van_der_waals_bond()
+        vdw_bonds = [b for b in mol.get_all_edges() if b.is_van_der_waals()]
+        self.assertEqual(len(vdw_bonds), 1)
+        vdw_bond = vdw_bonds[0]
+        self.assertEqual(vdw_bond.atom1.symbol,'N')
+
+
 ################################################################################
 
 if __name__ == '__main__':
diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd
index 586c4a6796..77ec42bfc2 100644
--- a/rmgpy/reaction.pxd
+++ b/rmgpy/reaction.pxd
@@ -58,6 +58,9 @@ cdef class Reaction:
     cdef public bint is_forward
     cdef public bint allow_max_rate_violation
     cdef public object rank
+    cdef public str metal
+    cdef public str facet
+    cdef public str site
 
     cpdef bint is_isomerization(self)
 
diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index dfb598784f..c0345f82d1 100644
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -53,7 +53,7 @@
 from rmgpy.exceptions import ReactionError, KineticsError
 from rmgpy.kinetics import KineticsData, ArrheniusBM, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, \
     PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \
-    SurfaceArrheniusBEP, StickingCoefficientBEP
+    SurfaceArrheniusBEP, StickingCoefficientBEP, SurfaceArrheniusBM
 from rmgpy.kinetics.arrhenius import Arrhenius  # Separate because we cimport from rmgpy.kinetics.arrhenius
 from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient  # Separate because we cimport from rmgpy.kinetics.surface
 from rmgpy.kinetics.diffusionLimited import diffusion_limiter
@@ -113,6 +113,9 @@ def __init__(self,
                  allow_max_rate_violation=False,
                  rank=None,
                  comment='',
+                 metal=None,
+                 facet=None,
+                 site=None,
                  is_forward=None,
                  ):
         self.index = index
@@ -134,6 +137,9 @@ def __init__(self,
         self.is_forward = is_forward
         self.allow_max_rate_violation = allow_max_rate_violation
         self.rank = rank
+        self.metal = metal
+        self.facet = facet
+        self.site = site
 
     def __repr__(self):
         """
@@ -157,7 +163,10 @@ def __repr__(self):
         if self.elementary_high_p: string += 'elementary_high_p={0}, '.format(self.elementary_high_p)
         if self.comment != '': string += 'comment={0!r}, '.format(self.comment)
         if self.rank is not None: string += 'rank={0!r},'.format(self.rank)
-        string = string[:-2] + ')'
+        if self.metal is not None: string += 'metal={0!r},'.format(self.metal)
+        if self.facet is not None: string += 'facet={0!r},'.format(self.facet)
+        if self.site is not None: string += 'site={0!r},'.format(self.site)
+        string = string[:-1] + ')'
         return string
 
     def __str__(self):
@@ -197,8 +206,12 @@ def __reduce__(self):
                            self.pairs,
                            self.allow_pdep_route,
                            self.elementary_high_p,
+                           self.allow_max_rate_violation,
                            self.rank,
-                           self.comment
+                           self.comment,
+                           self.metal,
+                           self.facet,
+                           self.site,
                            ))
 
     @property
@@ -231,6 +244,16 @@ def degeneracy(self, new):
         # set new degeneracy
         self._degeneracy = new
 
+    def get_metal_label(self):
+        """
+        retrieve the metal label (str) for the reaction
+        """
+        if self.facet is None:
+            metal = self.metal # could be None
+        else:
+            metal = self.metal + self.facet
+        return metal
+
     def to_chemkin(self, species_list=None, kinetics=True):
         """
         Return the chemkin-formatted string for this reaction.
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index 181024457e..8521a9b596 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -101,6 +101,7 @@ def database(
 def catalyst_properties(bindingEnergies=None,
                         surfaceSiteDensity=None,
                         metal=None,
+                        facet=None,
                         coverageDependence=False):
     """
     Specify the properties of the catalyst.
@@ -119,12 +120,19 @@ def catalyst_properties(bindingEnergies=None,
                          "or the surfaceSiteDensity and bindingEnergies, but not both.")
 
     if metal:
+        if facet:
+            metal_label = metal + facet
+        else:
+            metal_label = metal
         try:
-            logging.info("Using catalyst surface properties from metal %r.", metal)
-            rmg.binding_energies = metal_db.get_binding_energies(metal)
-            rmg.surface_site_density = metal_db.get_surface_site_density(metal)
+            logging.info("Using catalyst surface properties from metal %r.", metal_label)
+            rmg.metal = metal
+            rmg.facet = facet
+
+            rmg.binding_energies = metal_db.get_binding_energies(metal_label)
+            rmg.surface_site_density = metal_db.get_surface_site_density(metal_label)
         except DatabaseError:
-            logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal)
+            logging.error('Metal %r missing from surface library. Please specify both metal and facet.', metal_label)
             raise
     else: # metal not specified
         if bindingEnergies is None:
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index bd7f8f896d..8e5a706a49 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -184,6 +184,8 @@ def clear(self):
         self.kinetics_estimator = 'group additivity'
         self.solvent = None
         self.diffusion_limiter = None
+        self.metal = None
+        self.facet = None
         self.surface_site_density = None
         self.binding_energies = None
         self.coverage_dependence = False
@@ -268,7 +270,12 @@ def load_input(self, path=None):
             self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
             self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
         self.reaction_model.coverage_dependence = self.coverage_dependence
-            
+
+        if self.metal:
+            self.reaction_model.metal = self.metal
+        if self.facet:
+            self.reaction_model.facet = self.facet
+
         self.reaction_model.verbose_comments = self.verbose_comments
         self.reaction_model.save_edge_species = self.save_edge_species
 
@@ -382,6 +389,10 @@ def load_database(self):
             depository=False,  # Don't bother loading the depository information, as we don't use it
         )
 
+        # add metal attrs to kinetics database
+        self.database.kinetics.metal = self.metal
+        self.database.kinetics.facet = self.facet
+
         # Turn off reversibility for families with three products if desired
         if not self.trimolecular_product_reversible:
             for family in self.database.kinetics.families.values():
@@ -448,6 +459,7 @@ def load_database(self):
             for family in self.database.kinetics.families.values():
                 if not family.auto_generated:
                     family.fill_rules_by_averaging_up(verbose=self.verbose_comments)
+                family.split_template()
 
     def initialize(self, **kwargs):
         """
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 90d73ea989..a3bc0c319b 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -235,6 +235,8 @@ def __init__(self, core=None, edge=None, surface=None):
             3  R!H u1 px c0 {2,S}
             4  H u0 p0  c0 {1,S}
             """)]
+        self.metal = None
+        self.facet = None
 
     def check_for_existing_species(self, molecule):
         """
diff --git a/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py b/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py
index 720f9d1443..2560cbaa8a 100644
--- a/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py
+++ b/rmgpy/test_data/testing_database/kinetics/families/intra_H_migration/groups.py
@@ -11,6 +11,10 @@
 
 reversible = True
 
+reactantNum = 1
+
+productNum = 1
+
 autoGenerated = False
 
 recipe(actions=[
diff --git a/testing/databaseTest.py b/testing/databaseTest.py
index 8be0a6e7eb..82d8521ecb 100644
--- a/testing/databaseTest.py
+++ b/testing/databaseTest.py
@@ -242,12 +242,12 @@ def test_thermo(self):
             yield test, group_name
 
             # tests for adsorption groups
-            if 'adsorption' in group_name.lower():
-                test = lambda x: self.check_surface_thermo_groups_have_surface_attributes(group_name, group)
-                test_name = "Thermo surface groups {0}: Entry has metal attributes?".format(group_name)
-                test.description = test_name
-                self.compat_func_name = test_name
-                yield test, group_name
+            # if 'adsorption' in group_name.lower():
+            #     test = lambda x: self.check_surface_thermo_groups_have_surface_attributes(group_name, group)
+            #     test_name = "Thermo surface groups {0}: Entry has metal attributes?".format(group_name)
+            #     test.description = test_name
+            #     self.compat_func_name = test_name
+            #     yield test, group_name
 
         for library_name, library in self.database.thermo.libraries.items():
             if 'surface' in library_name.lower():
@@ -1514,30 +1514,30 @@ def make_error_message(reactants, message=''):
             boo = True
         if boo:
             raise ValueError("Error Occurred. See log for details.")
-    def check_surface_thermo_groups_have_surface_attributes(self, group_name, group):
-        """
-        Tests that each entry in the surface thermo groups has a 'metal' and 'facet' attribute, 
-        describing which metal the data came from.
-        """
-        failed = False
-        for entry in group.entries.values():
-            if isinstance(entry.data, rmgpy.thermo.thermodata.ThermoData):
-                if 'Pt' in group_name:
-                    if entry.metal is not 'Pt':
-                        logging.error(f'Expected {entry} metal attribute in {group_name} group to match Pt, but was {entry.metal}')
-                        failed = True
-                if '111' in group_name:
-                    if entry.facet is not '111':
-                        logging.error(f'Expected {entry} facet attribute in {group_name} group to match 111, but was {entry.facet}')
-                        failed = True
-                if not entry.metal:
-                    logging.error(f'Expected a metal attribute for {entry} in {group_name} group but found {entry.metal!r}')
-                    failed = True
-                if not entry.facet:
-                    logging.error(f'Expected a facet attribute for {entry} in {group_name} group but found {entry.facet!r}')
-                    failed = True
-        if failed:
-            raise ValueError("Error occured in databaseTest. Please check log warnings for all error messages.")
+    # def check_surface_thermo_groups_have_surface_attributes(self, group_name, group):
+    #     """
+    #     Tests that each entry in the surface thermo groups has a 'metal' and 'facet' attribute, 
+    #     describing which metal the data came from.
+    #     """
+    #     failed = False
+    #     for entry in group.entries.values():
+    #         if isinstance(entry.data, rmgpy.thermo.thermodata.ThermoData):
+    #             if 'Pt' in group_name:
+    #                 if entry.metal is not 'Pt':
+    #                     logging.error(f'Expected {entry} metal attribute in {group_name} group to match Pt, but was {entry.metal}')
+    #                     failed = True
+    #             if '111' in group_name:
+    #                 if entry.facet is not '111':
+    #                     logging.error(f'Expected {entry} facet attribute in {group_name} group to match 111, but was {entry.facet}')
+    #                     failed = True
+    #             if not entry.metal:
+    #                 logging.error(f'Expected a metal attribute for {entry} in {group_name} group but found {entry.metal!r}')
+    #                 failed = True
+    #             if not entry.facet:
+    #                 logging.error(f'Expected a facet attribute for {entry} in {group_name} group but found {entry.facet!r}')
+    #                 failed = True
+    #     if failed:
+    #         raise ValueError("Error occured in databaseTest. Please check log warnings for all error messages.")
 
     def check_surface_thermo_libraries_have_surface_attributes(self, library_name, library):
         """