diff --git a/documentation/source/reference/molecule/adjlist.rst b/documentation/source/reference/molecule/adjlist.rst
index a259682a11..7a4980ae07 100644
--- a/documentation/source/reference/molecule/adjlist.rst
+++ b/documentation/source/reference/molecule/adjlist.rst
@@ -34,8 +34,8 @@ alphanumeric characters and the underscore, as if an identifier in many popular
 programming languages. However, strictly speaking any non-space ASCII character
 is allowed.
 
-The subsequent lines may contain keyword-value pairs. Currently there is only
-one keyword, ``multiplicity``.
+The subsequent lines may contain keyword-value pairs. Currently there can be: 
+``multiplicity``, ``metal`` and ``facet``.
 
 For species or molecule declarations, the value after ``multiplicity`` defines 
 the spin multiplicity of the molecule. E.g. ``multiplicity 1`` for most ground state 
@@ -44,7 +44,7 @@ and ``multiplicity 3`` for a triplet biradical.
 If the ``multiplicity`` line is not present then a value of 
 (1 + number of unpaired electrons) is assumed. 
 Thus, it can usually be omitted, but if present can be used to distinguish,
-for example, singlet CH2 from triplet CH2.
+for example, singlet CH2 from triplet CH2. 
 
 If defining a Functional :class:`~rmgpy.molecule.Group`, then the value must be a list,
 which defines the multiplicities that will be matched by the group, eg. 
@@ -53,6 +53,9 @@ If a wildcard is desired, the line ``'multiplicity x`` can be used instead to ac
 all multiplicities.  If the multiplicity line is omitted altogether, then a wildcard 
 is assumed.
 
+``metal`` and ``facet`` work similarly and will correspond to lines like ``metal Fe``, 
+``metal [Fe,Cu,Ag]``, ``facet 111``, ``facet [111,211,110]``. 
+
 e.g. the following two group adjlists represent identical groups. ::
 
     group1
@@ -69,7 +72,7 @@ each subsequent line describes a single atom and its
 local bond structure. The format of these lines is a whitespace-delimited list
 with tokens ::
 
-    <number> [<label>] <element> u<unpaired> [p<pairs>] [c<charge>] <bondlist>
+    <number> [<label>] <element> u<unpaired> [p<pairs>] [c<charge>] [s<site>] [m<morphology>] <bondlist>
 
 The first item is the number used to identify that atom. Any number may be used,
 though it is recommended to number the atoms sequentially starting from one.
@@ -86,6 +89,8 @@ a sequence of tokens describing the electronic state of the atom:
 * ``p0`` number of lone **pairs** of electrons, common on oxygen and nitrogen.
 * ``c0`` formal **charge** on the atom, e.g. ``c-1`` (negatively charged),
   ``c0``, ``c+1`` (positively charged)
+* ``s`` the site type a site atom is e.g. ``s"fcc"``
+* ``m`` the morphology of a site atom e.g. ``m"terrace"``
 
 For :class:`~rmgpy.molecule.Molecule` definitions:
 The value must be a single integer (and for charge must have a + or - sign if not equal to 0)
@@ -94,7 +99,7 @@ The number of lone pairs and the formal charge are assumed to be zero if omitted
 
 For :class:`~rmgpy.molecule.Group` definitions:
 The value can be an integer or a list of integers (with signs, for charges),
-eg. ``u[0,1,2]`` or ``c[0,+1,+2,+3,+4]``, or may be a wildcard ``x`` 
+eg. ``u[0,1,2]``, ``c[0,+1,+2,+3,+4]``, or ``s["hcp","fcc"] or may be a wildcard ``x`` 
 which matches any valid value, 
 eg. ``px`` is the same as ``p[0,1,2,3,4, ...]`` and ``cx`` is the same as
 ``c[...,-4,-3,-2,-1,0,+1,+2,+3,+4,...]``. Lists must be enclosed is square brackets,
@@ -169,6 +174,12 @@ The allowed element types, radicals, and bonds are listed in the following table
  |                      | T        | Triple bond         |
  |                      +----------+---------------------+
  |                      | B        | Benzene bond        |
+ |                      +--------------------------------+
+ |                      | vdW      | Van der Waals bond  |
+ |                      +--------------------------------+
+ |                      | H        | Hydrogen bond       |
+ |                      +--------------------------------+
+ |                      | R        | Reaction bond       |
  +----------------------+----------+---------------------+
 
 
diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py
index 92270b6595..f559dda933 100644
--- a/rmgpy/molecule/adjlist.py
+++ b/rmgpy/molecule/adjlist.py
@@ -438,8 +438,10 @@ def from_old_adjacency_list(adjlist, group=False, saturate_h=False):
     except InvalidAdjacencyListError:
         logging.error("Troublesome adjacency list:\n" + adjlist)
         raise
-
-    return atoms, multiplicity
+    if group:
+        return atoms, multiplicity, [], []
+    else:
+        return atoms, multiplicity, '', ''
 
 
 ###############################
@@ -460,7 +462,7 @@ def from_old_adjacency_list(adjlist, group=False, saturate_h=False):
                             r'\s*$')  # the end!
 
 
-def from_adjacency_list(adjlist, group=False, saturate_h=False):
+def from_adjacency_list(adjlist, group=False, saturate_h=False, check_consistency=True):
     """
     Convert a string adjacency list `adjlist` into a set of :class:`Atom` and
     :class:`Bond` objects.
@@ -524,11 +526,17 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
             multiplicity = int(line.split()[1])
         if len(lines) == 0:
             raise InvalidAdjacencyListError('No atoms specified in adjacency list: \n{0}'.format(adjlist))
-
+        
     mistake1 = re.compile(r'\{[^}]*\s+[^}]*\}')
+    if group:
+        metal = []
+        facet = []
+    else:
+        metal = ''
+        facet = ''
     # Iterate over the remaining lines, generating Atom or GroupAtom objects
     for line in lines:
-
+            
         # Sometimes people put spaces after commas, which messes up the
         # parse-by-whitespace. Examples include '[Cd, Ct]'.
         if mistake1.search(line):
@@ -545,7 +553,49 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
         # Skip if blank line
         if len(data) == 0:
             continue
-
+        
+        if line.split()[0] == 'metal':
+            if group:
+                match = re.match(r'\s*metal\s+\[\s*[\w,\s]+\s*\]\s*$', line)
+                if not match:
+                    rematch = re.match(r'\s*metal\s+x\s*$', line)
+                    if not rematch:
+                        raise InvalidAdjacencyListError("Invalid metal line '{0}'. Should be a list like "
+                                                        "'metal [Cu,Fe,Ag]' or a wildcard 'metal x'".format(line))
+                else:
+                    out = line.split('[')[1][:-1]
+                    metal = [x.strip() for x in out.split(',') if x.strip() != '']
+        
+            else:
+                match = re.match(r'\s*metal\s+\w+\s*$', line)
+                if not match:
+                    raise InvalidAdjacencyListError("Invalid metal line '{0}'. Should be a string like "
+                                                "'metal Fe'".format(line))
+                metal = line.split()[1].strip()
+            
+            continue
+                
+        if line.split()[0] == 'facet':
+            if group:
+                match = re.match(r'\s*facet\s+\[\s*[\w,\s]+\s*\]\s*$', line)
+                if not match:
+                    rematch = re.match(r'\s*facet\s+x\s*$', line)
+                    if not rematch:
+                        raise InvalidAdjacencyListError("Invalid facet line '{0}'. Should be a list like "
+                                                        "'facet [111,211,110]' or a wildcard 'facet x'".format(line))
+                else:
+                    out = line.split('[')[1][:-1]
+                    facet = [x.strip() for x in out.split(',') if x.strip() != '']
+        
+            else:
+                match = re.match(r'\s*facet\s+\w+\s*$', line)
+                if not match:
+                    raise InvalidAdjacencyListError("Invalid facet line '{0}'. Should be a string like "
+                                                "'facet 111'".format(line))
+                facet = line.split()[1].strip()
+            
+            continue
+        
         # First item is index for atom
         # Sometimes these have a trailing period (as if in a numbered list),
         # so remove it just in case
@@ -677,6 +727,32 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
             if not group:
                 partial_charges.append(0)
 
+        # Next the sites (if provided)
+        sites = []
+        if len(data) > index:
+            s_state = data[index]
+            if s_state[0] == 's':
+                if s_state[1] == '[':
+                    s_state = s_state[2:-1].split(',')
+                else:
+                    s_state = [s_state[1:]]
+                for s in s_state:
+                    sites.append(s[1:-1])
+                index += 1
+            
+        # Next the morphologys (if provided)
+        morphologies = []
+        if len(data) > index:
+            m_state = data[index]
+            if m_state[0] == 'm':
+                if m_state[1] == '[':
+                    m_state = m_state[2:-1].split(',')
+                else:
+                    m_state = [m_state[1:]]
+                for m in m_state:
+                    morphologies.append(m[1:-1])
+                index += 1
+        
         # Next the isotope (if provided)
         isotope = -1
         if len(data) > index:
@@ -695,12 +771,20 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
 
         # Create a new atom based on the above information
         if group:
-            atom = GroupAtom(atom_type, unpaired_electrons, partial_charges, label, lone_pairs, props)
+            atom = GroupAtom(atom_type, unpaired_electrons, partial_charges, label, lone_pairs, sites, morphologies, props)
         else:
             # detect if this is cutting label or atom
             _ , cutting_label_list = Fragment().detect_cutting_label(atom_type[0])
             if cutting_label_list == []:
-                atom = Atom(atom_type[0], unpaired_electrons[0], partial_charges[0], label, lone_pairs[0])
+                if sites == []:
+                    site = ''
+                else:
+                    site = sites[0]
+                if morphologies == []:
+                    morphology = ''
+                else:
+                    morphology = morphologies[0]
+                atom = Atom(atom_type[0], unpaired_electrons[0], partial_charges[0], label, lone_pairs[0], site, morphology)
                 if isotope != -1:
                     atom.element = get_element(atom.number, isotope)
             else:
@@ -772,7 +856,7 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
             Saturator.saturate(atoms)
 
     # Consistency checks
-    if not group:
+    if not group and check_consistency:
         # Molecule consistency check
         # Electron and valency consistency check for each atom
         for atom in atoms:
@@ -787,13 +871,19 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False):
         ConsistencyChecker.check_multiplicity(n_rad, multiplicity)
         for atom in atoms:
             ConsistencyChecker.check_hund_rule(atom, multiplicity)
-        return atoms, multiplicity
+        return atoms, multiplicity, metal, facet
     else:
         # Currently no group consistency check
-        return atoms, multiplicity
+        if not group:
+            if multiplicity is None:
+                n_rad = sum([atom.radical_electrons for atom in atoms])
+                multiplicity = n_rad + 1
+
+        return atoms, multiplicity, metal, facet
+
 
 
-def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=False, remove_lone_pairs=False,
+def to_adjacency_list(atoms, multiplicity, metal='', facet='', label=None, group=False, remove_h=False, remove_lone_pairs=False,
                       old_style=False):
     """
     Convert a chemical graph defined by a list of `atoms` into a string
@@ -821,10 +911,18 @@ def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=Fal
             # Functional group should have a list of possible multiplicities.  
             # If the list is empty, then it does not need to be written
             adjlist += 'multiplicity [{0!s}]\n'.format(','.join(str(i) for i in multiplicity))
+        if metal:
+            adjlist += 'metal [{0!s}]\n'.format(','.join(i for i in metal))
+        if facet:
+            adjlist += 'facet [{0!s}]\n'.format(','.join(i for i in facet))
     else:
         assert isinstance(multiplicity, int), "Molecule should have an integer multiplicity"
         if multiplicity != 1 or any(atom.radical_electrons for atom in atoms):
             adjlist += 'multiplicity {0!r}\n'.format(multiplicity)
+        if metal:
+            adjlist += f"metal {metal}\n"
+        if facet:
+            adjlist += f"facet {facet}\n"
 
     # Determine the numbers to use for each atom
     atom_numbers = {}
@@ -843,6 +941,8 @@ def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=Fal
     atom_charge = {}
     atom_isotope = {}
     atom_props = {}
+    atom_site = {}
+    atom_morphology = {}
     if group:
         for atom in atom_numbers:
             # Atom type(s)
@@ -874,6 +974,22 @@ def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=Fal
             else:
                 atom_charge[atom] = '[{0}]'.format(','.join(['+'+str(charge) if charge > 0 else ''+str(charge) for charge in atom.charge]))
 
+            # Sites
+            if len(atom.site) == 1:
+                atom_site[atom] = "\"" + atom.site[0] + "\""
+            elif len(atom.site) == 0:
+                atom_site[atom] = None  # Empty list indicates wildcard
+            else:
+                atom_site[atom] = '["{0}"]'.format('","'.join(s for s in atom.site))
+            
+            # Morphologies
+            if len(atom.morphology) == 1:
+                atom_morphology[atom] = "\"" + atom.morphology[0] + "\""
+            elif len(atom.morphology) == 0:
+                atom_morphology[atom] = None  # Empty list indicates wildcard
+            else:
+                atom_morphology[atom] = '["{0}"]'.format('","'.join(s for s in atom.morphology))
+                
             # Isotopes
             atom_isotope[atom] = -1
 
@@ -892,6 +1008,16 @@ def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=Fal
             atom_lone_pairs[atom] = str(atom.lone_pairs)
             # Partial Charge(s)
             atom_charge[atom] = '+' + str(atom.charge) if atom.charge > 0 else '' + str(atom.charge)
+            # Sites
+            if atom.site:
+                atom_site[atom] = "\"" + atom.site + "\""
+            else:
+                atom_site[atom] = None
+            # Morphology
+            if atom.morphology:
+                atom_morphology[atom] = "\"" + atom.morphology + "\""
+            else:
+                atom_morphology[atom] = None
             # Isotopes
             if isinstance(atom, Atom):
                 atom_isotope[atom] = atom.element.isotope
@@ -926,6 +1052,12 @@ def to_adjacency_list(atoms, multiplicity, label=None, group=False, remove_h=Fal
         # Partial charges
         if atom_charge[atom] is not None:
             adjlist += ' c{0}'.format(atom_charge[atom])
+        # Sites
+        if atom_site[atom]:
+            adjlist += ' s{0}'.format(atom_site[atom])
+        # Morphologies
+        if atom_morphology[atom]:
+            adjlist += ' m{0}'.format(atom_morphology[atom])
         # Isotopes
         if atom_isotope[atom] != -1:
             adjlist += ' i{0}'.format(atom_isotope[atom])
diff --git a/rmgpy/molecule/adjlistTest.py b/rmgpy/molecule/adjlistTest.py
index ad6585e8b5..3366f4f0d8 100644
--- a/rmgpy/molecule/adjlistTest.py
+++ b/rmgpy/molecule/adjlistTest.py
@@ -171,6 +171,58 @@ def test_atom_props(self):
         adjlist2 = group.to_adjacency_list()
 
         self.assertEqual(adjlist.strip(), adjlist2.strip())
+        
+    def test_metal_facet_site_morphology(self):
+        adjlist1 = """metal Cu
+facet 111
+1 *3 X u0 p0 c0 s"atop" m"terrace"  {2,S} {4,S}
+2 *1 O u0 p2 c0 {1,S} {3,R}
+3 *2 H u0 p0 c0 {2,R} {4,R}
+4 *4 X u0 p0 c0 s"hcp" m"terrace" {3,R} {1,S}"""
+        
+        adjlist2 = """multiplicity [1]
+metal [Cu, Fe, CuO2 ]
+facet [111, 211, 1101, 110, ]
+1 *3 X u0 p0 c0 s["atop","fcc"] m"terrace"  {2,S} {4,S}
+2 *1 O u0 p2 c0 {1,S} {3,R}
+3 *2 H u0 p0 c0 {2,R} {4,R}
+4 *4 X u0 p0 c0 s"hcp" m["terrace","sc"] {3,R} {1,S}"""
+        
+        adjlist1test = """metal Cu
+facet 111
+1 *3 X u0 p0 c0 s"atop" m"terrace" {2,S} {4,S}
+2 *1 O u0 p2 c0 {1,S} {3,R}
+3 *2 H u0 p0 c0 {2,R} {4,R}
+4 *4 X u0 p0 c0 s"hcp" m"terrace" {1,S} {3,R}"""
+        
+        adjlist2test = """multiplicity [1]
+metal [Cu,Fe,CuO2]
+facet [111,211,1101,110]
+1 *3 X u0 p0 c0 s["atop","fcc"] m"terrace" {2,S} {4,S}
+2 *1 O u0 p2 c0 {1,S} {3,R}
+3 *2 H u0 p0 c0 {2,R} {4,R}
+4 *4 X u0 p0 c0 s"hcp" m["terrace","sc"] {1,S} {3,R}"""
+        mol = Molecule().from_adjacency_list(adjlist1,check_consistency=False)
+        group = Group().from_adjacency_list(adjlist2)
+        
+        self.assertEqual(mol.metal,"Cu")
+        self.assertEqual(mol.facet,"111")
+        self.assertEqual(group.metal,["Cu","Fe","CuO2"])
+        self.assertEqual(group.facet,["111","211","1101","110"])
+        
+        self.assertEqual(mol.atoms[0].site,"atop")
+        self.assertEqual(mol.atoms[3].site,"hcp")
+        self.assertEqual(mol.atoms[0].morphology, "terrace")
+        self.assertEqual(mol.atoms[3].morphology, "terrace")
+        
+        self.assertEqual(group.atoms[0].site,["atop","fcc"])
+        self.assertEqual(group.atoms[3].site,["hcp"])
+        self.assertEqual(group.atoms[0].morphology, ["terrace"])
+        self.assertEqual(group.atoms[3].morphology,["terrace","sc"])
+        
+        self.assertEqual(mol.to_adjacency_list().strip(),adjlist1test.strip())
+        self.assertEqual(group.to_adjacency_list().strip(), adjlist2test.strip())
+
 
 
 class TestMoleculeAdjLists(unittest.TestCase):
diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py
index 41620ee322..6b32247e45 100644
--- a/rmgpy/molecule/atomtype.py
+++ b/rmgpy/molecule/atomtype.py
@@ -249,8 +249,8 @@ def get_features(self):
                              benzene=[0], lone_pairs=[0])
 # Occupied surface site:
 ATOMTYPES['Xo']   = AtomType('Xo', generic=['X'], specific=[],
-                             single=[0, 1], all_double=[0, 1], r_double=[], o_double=[], s_double=[], triple=[0, 1],
-                             quadruple=[0, 1], benzene=[0], lone_pairs=[0])
+                             single=[0, 1, 2, 3, 4, 5, 6, 7, 8], all_double=[0, 1, 2, 3, 4, 5, 6, 7, 8], r_double=[], o_double=[], s_double=[], triple=[0, 1, 2, 3, 4, 5, 6, 7, 8],
+                             quadruple=[0, 1, 2, 3, 4, 5, 6, 7, 8], benzene=[0], lone_pairs=[0])
 
 # Non-surface atomTypes, R being the most generic:
 ATOMTYPES['R']    = AtomType(label='R', generic=[], specific=[
diff --git a/rmgpy/molecule/fragment.py b/rmgpy/molecule/fragment.py
index 298b1ae24b..f7742e67dd 100644
--- a/rmgpy/molecule/fragment.py
+++ b/rmgpy/molecule/fragment.py
@@ -28,6 +28,8 @@ def __init__(self, name='', label='', id=-1):
         self.isotope = -1
         self.id = id
         self.mass = 0
+        self.site = ''
+        self.morphology = ''
 
     def __str__(self):
         """
@@ -148,6 +150,8 @@ def __init__(self,
         self.props = props or {}
         self.multiplicity = multiplicity
         self.reactive = reactive
+        self.metal = ''
+        self.facet = ''
 
         if inchi and smiles:
             logging.warning('Both InChI and SMILES provided for Fragment instantiation, '
@@ -1290,7 +1294,7 @@ def from_adjacency_list(self, adjlist, saturate_h=False, raise_atomtype_exceptio
         """
         from rmgpy.molecule.adjlist import from_adjacency_list
         
-        self.vertices, self.multiplicity = from_adjacency_list(adjlist, group=False, saturate_h=saturate_h)
+        self.vertices, self.multiplicity, self.site, self.morphology = from_adjacency_list(adjlist, group=False, saturate_h=saturate_h)
         self.update_atomtypes(raise_exception=raise_atomtype_exception)
         
         # Check if multiplicity is possible
diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd
index fc0efbdf38..e74f6bea7d 100644
--- a/rmgpy/molecule/group.pxd
+++ b/rmgpy/molecule/group.pxd
@@ -38,13 +38,16 @@ cdef class GroupAtom(Vertex):
     cdef public list charge
     cdef public str label
     cdef public list lone_pairs
-
+    cdef public list site 
+    cdef public list morphology
     cdef public dict props
 
     cdef public list reg_dim_atm
     cdef public list reg_dim_u
     cdef public list reg_dim_r
-
+    cdef public list reg_dim_site
+    cdef public list reg_dim_morphology
+    
     cpdef Vertex copy(self)
 
     cpdef _change_bond(self, short order)
@@ -126,7 +129,8 @@ cdef class Group(Graph):
 
     cdef public dict props
     cdef public list multiplicity
-
+    cdef public list metal 
+    cdef public list facet
     # These read-only attribues act as a "fingerprint" for accelerating
     # subgraph isomorphism checks
     cdef public dict elementCount
@@ -166,7 +170,7 @@ cdef class Group(Graph):
 
     cpdef dict get_element_count(self)
 
-    cpdef from_adjacency_list(self, str adjlist)
+    cpdef from_adjacency_list(self, str adjlist, bint check_consistency=?)
 
     cpdef to_adjacency_list(self, str label=?)
     
diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py
index 20b95c3953..1866b9b644 100644
--- a/rmgpy/molecule/group.py
+++ b/rmgpy/molecule/group.py
@@ -63,10 +63,14 @@ class GroupAtom(Vertex):
     `label`              ``str``             A string label that can be used to tag individual atoms
     `lone_pairs`         ``list``            The number of lone electron pairs
     `charge`             ``list``            The partial charge of the atom
+    `site`               ``list``            The allowed adsorption sites
+    `morphology`         ``list``            The allowed morphologies
     `props`              ``dict``            Dictionary for storing additional atom properties
     `reg_dim_atm`        ``list``            List of atom types that are free dimensions in tree optimization
     `reg_dim_u`          ``list``            List of unpaired electron numbers that are free dimensions in tree optimization
     `reg_dim_r`          ``list``            List of inRing values that are free dimensions in tree optimization
+    `reg_dim_site`       ``list``            List of sites that are free dimensions in tree optimization
+    `reg_dim_morphology` ``list``            List of morphologies that are free dimensions in tree optimization
     ==================== =================== ====================================
 
     Each list represents a logical OR construct, i.e. an atom will match the
@@ -76,7 +80,8 @@ class GroupAtom(Vertex):
     order to match.
     """
 
-    def __init__(self, atomtype=None, radical_electrons=None, charge=None, label='', lone_pairs=None, props=None):
+    def __init__(self, atomtype=None, radical_electrons=None, charge=None, label='', lone_pairs=None, site=None, morphology=None, 
+                 props=None):
         Vertex.__init__(self)
         self.atomtype = atomtype or []
         for index in range(len(self.atomtype)):
@@ -86,12 +91,15 @@ def __init__(self, atomtype=None, radical_electrons=None, charge=None, label='',
         self.charge = charge or []
         self.label = label
         self.lone_pairs = lone_pairs or []
-
+        self.site = site or []
+        self.morphology = morphology or []
         self.props = props or {}
 
         self.reg_dim_atm = [[], []]
         self.reg_dim_u = [[], []]
         self.reg_dim_r = [[], []]
+        self.reg_dim_site = [[], []]
+        self.reg_dim_morphology = [[], []]
 
     def __reduce__(self):
         """
@@ -107,7 +115,8 @@ def __reduce__(self):
         atomtype = self.atomtype
         if atomtype is not None:
             atomtype = [a.label for a in atomtype]
-        return (GroupAtom, (atomtype, self.radical_electrons, self.charge, self.label, self.lone_pairs, self.props), d)
+        return (GroupAtom, (atomtype, self.radical_electrons, self.charge, self.label, self.lone_pairs, self.site, 
+                            self.morphology, self.props), d)
 
     def __setstate__(self, d):
         """
@@ -146,6 +155,8 @@ def copy(self):
             self.charge[:],
             self.label,
             self.lone_pairs[:],
+            self.site[:],
+            self.morphology[:],
             deepcopy(self.props),
         )
 
@@ -404,6 +415,32 @@ def equivalent(self, other, strict=True):
                     if charge1 == charge2: break
                 else:
                     return False
+        # Each site in self must have an equivalent in other (and vice versa)
+        for site1 in self.site:
+            if group.site:
+                for site2 in group.site:
+                    if site1 == site2: break
+                else:
+                    return False
+        for site1 in group.site:
+            if self.site:
+                for site2 in self.site:
+                    if site1 == site2: break
+                else:
+                    return False
+        # Each morphology in self must have an equivalent in other (and vice versa)
+        for morphology1 in self.morphology:
+            if group.morphology:
+                for morphology2 in group.morphology:
+                    if morphology1 == morphology2: break
+                else:
+                    return False
+        for morphology1 in group.morphology:
+            if self.morphology:
+                for morphology2 in self.morphology:
+                    if morphology1 == morphology2: break
+                else:
+                    return False
         # Other properties must have an equivalent in other (and vice versa)
         # Absence of the 'inRing' prop indicates a wildcard
         if 'inRing' in self.props and 'inRing' in group.props:
@@ -427,7 +464,8 @@ def is_specific_case_of(self, other):
         group = other
 
         cython.declare(atomType1=AtomType, atomtype2=AtomType, radical1=cython.short, radical2=cython.short,
-                       lp1=cython.short, lp2=cython.short, charge1=cython.short, charge2=cython.short)
+                       lp1=cython.short, lp2=cython.short, charge1=cython.short, charge2=cython.short,
+                       site1=str, site2=str, morphology1=str, morphology2=str)
         # Compare two atom groups for equivalence
         # Each atom type in self must have an equivalent in other (and vice versa)
         for atomType1 in self.atomtype:  # all these must match
@@ -464,6 +502,26 @@ def is_specific_case_of(self, other):
                         return False
         else:
             if group.charge: return False
+        # Each site in self must have an equivalent in other
+        if self.site:
+            for site1 in self.site:
+                if group.site:
+                    for site2 in group.site:
+                        if site1 == site2: break
+                    else:
+                        return False
+        else:
+            if group.site: return False
+        # Each morphology in self must have an equivalent in other
+        if self.morphology:
+            for morphology1 in self.morphology:
+                if group.morphology:
+                    for morphology2 in group.morphology:
+                        if morphology1 == morphology2: break
+                    else:
+                        return False
+        else:
+            if group.morphology: return False
         # Other properties must have an equivalent in other
         # Absence of the 'inRing' prop indicates a wildcard
         if 'inRing' in self.props and 'inRing' in group.props:
@@ -750,6 +808,8 @@ def get_order_str(self):
                 values.append('vdW')
             elif value == 0.1:
                 values.append('H')
+            elif value == 0.05:
+                values.append('R')
             else:
                 raise TypeError('Bond order number {} is not hardcoded as a string'.format(value))
         return values
@@ -775,6 +835,8 @@ def set_order_str(self, new_order):
                 values.append(1.5)
             elif value == 'H':
                 values.append(0.1)
+            elif value == 'R':
+                values.append(0.05)
             else:
                 # try to see if an float disguised as a string was input by mistake
                 try:
@@ -903,12 +965,28 @@ def is_hydrogen_bond(self, wildcards=False):
         """
         if wildcards:
             for order in self.order:
-                if abs(order) <= 1e-9:
+                if abs(order - 0.1) <= 1e-9:
                     return True
             else:
                 return False
         else:
-            return abs(self.order[0]) <= 1e-9 and len(self.order) == 1
+            return abs(self.order[0] - 0.1) <= 1e-9 and len(self.order) == 1
+        
+    def is_reaction_bond(self, wildcards=False):
+        """
+        Return ``True`` if the bond represents a reaction bond or ``False`` if
+        not. If `wildcards` is ``False`` we return False anytime there is more
+        than one bond order, otherwise we return ``True`` if any of the options
+        are reaction bonds.
+        """
+        if wildcards:
+            for order in self.order:
+                if abs(order - 0.05) <= 1e-9:
+                    return True
+            else:
+                return False
+        else:
+            return abs(self.order[0] - 0.05) <= 1e-9 and len(self.order) == 1
 
     def _change_bond(self, order):
         """
@@ -1026,15 +1104,19 @@ class Group(Graph):
     `atoms`             ``list``            Aliases for the `vertices` storing :class:`GroupAtom`
     `multiplicity`      ``list``            Range of multiplicities accepted for the group
     `props`             ``dict``            Dictionary of arbitrary properties/flags classifying state of Group object 
+    `metal`             ``list``            List of metals accepted for the group
+    `facet`             ``list``            List of facets accepted for the group
     =================== =================== ====================================
 
     Corresponding alias methods to Molecule have also been provided.
     """
 
-    def __init__(self, atoms=None, props=None, multiplicity=None):
+    def __init__(self, atoms=None, props=None, multiplicity=None, metal=None, facet=None):
         Graph.__init__(self, atoms)
         self.props = props or {}
         self.multiplicity = multiplicity or []
+        self.metal = metal or []
+        self.facet = facet or []
         self.elementCount = {}
         self.radicalCount = -1
         self.update()
@@ -1744,14 +1826,14 @@ def get_element_count(self):
 
         return element_count
 
-    def from_adjacency_list(self, adjlist):
+    def from_adjacency_list(self, adjlist, check_consistency=True):
         """
         Convert a string adjacency list `adjlist` to a molecular structure.
         Skips the first line (assuming it's a label) unless `withLabel` is
         ``False``.
         """
         from rmgpy.molecule.adjlist import from_adjacency_list
-        self.vertices, multiplicity = from_adjacency_list(adjlist, group=True)
+        self.vertices, multiplicity, self.metal, self.facet = from_adjacency_list(adjlist, group=True, check_consistency=check_consistency)
         if multiplicity is not None:
             self.multiplicity = multiplicity
         self.update()
@@ -1762,7 +1844,7 @@ def to_adjacency_list(self, label=''):
         Convert the molecular structure to a string adjacency list.
         """
         from rmgpy.molecule.adjlist import to_adjacency_list
-        return to_adjacency_list(self.vertices, multiplicity=self.multiplicity, label=label, group=True)
+        return to_adjacency_list(self.vertices, multiplicity=self.multiplicity, metal=self.metal, facet=self.facet, label=label, group=True)
 
     def update_fingerprint(self):
         """
@@ -1825,7 +1907,7 @@ def is_subgraph_isomorphic(self, other, initial_map=None, generate_initial_map=F
         be a :class:`Group` object, or a :class:`TypeError` is raised.
         """
         cython.declare(group=Group)
-        cython.declare(mult1=cython.short, mult2=cython.short)
+        cython.declare(mult1=cython.short, mult2=cython.short, m1=str, m2=str)
         cython.declare(a=GroupAtom, L=list)
         # It only makes sense to compare a Group to a Group for subgraph
         # isomorphism, so raise an exception if this is not what was requested
@@ -1874,6 +1956,24 @@ def is_subgraph_isomorphic(self, other, initial_map=None, generate_initial_map=F
                         return False
         else:
             if group.multiplicity: return False
+        if self.metal:
+            for m1 in self.metal:
+                if group.metal:
+                    for m2 in group.metal:
+                        if m1 == m2: break
+                    else:
+                        return False
+        else:
+            if group.metal: return False
+        if self.facet:
+            for m1 in self.facet:
+                if group.facet:
+                    for m2 in group.facet:
+                        if m1 == m2: break
+                    else:
+                        return False
+        else:
+            if group.facet: return False
         # Do the isomorphism comparison
         return Graph.is_subgraph_isomorphic(self, other, initial_map, save_order=save_order)
 
@@ -1890,7 +1990,7 @@ def find_subgraph_isomorphisms(self, other, initial_map=None, save_order=False):
         :class:`TypeError` is raised.
         """
         cython.declare(group=Group)
-        cython.declare(mult1=cython.short, mult2=cython.short)
+        cython.declare(mult1=cython.short, mult2=cython.short, m1=str, m2=str)
 
         # It only makes sense to compare a Group to a Group for subgraph
         # isomorphism, so raise an exception if this is not what was requested
@@ -1909,7 +2009,27 @@ def find_subgraph_isomorphisms(self, other, initial_map=None, save_order=False):
         else:
             if group.multiplicity:
                 return []
-
+        if self.metal:
+            for m1 in self.metal:
+                if group.metal:
+                    for m2 in group.metal:
+                        if m1 == m2: break
+                    else:
+                        return []
+        else:
+            if group.metal:
+                return []
+        if self.facet:
+            for m1 in self.facet:
+                if group.facet:
+                    for m2 in group.facet:
+                        if m1 == m2: break
+                    else:
+                        return []
+        else:
+            if group.facet:
+                return []
+            
         # Do the isomorphism comparison
         return Graph.find_subgraph_isomorphisms(self, other, initial_map, save_order=save_order)
 
diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index f3cb1c8411..3d5062f8bb 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -45,6 +45,8 @@ cdef class Atom(Vertex):
     cdef public AtomType atomtype
     cdef public np.ndarray coords
     cdef public short lone_pairs
+    cdef public str site 
+    cdef public str morphology
     cdef public int id
     cdef public dict props
     
@@ -146,6 +148,8 @@ cdef class Molecule(Graph):
     cdef public int multiplicity
     cdef public bint reactive
     cdef public dict props
+    cdef public str metal
+    cdef public str facet
     cdef str _fingerprint
     cdef str _inchi
     cdef str _smiles
@@ -217,7 +221,7 @@ cdef class Molecule(Graph):
     cpdef from_smiles(self, str smilesstr, backend=?, bint raise_atomtype_exception=?)
 
     cpdef from_adjacency_list(self, str adjlist, bint saturate_h=?, bint raise_atomtype_exception=?,
-                              bint raise_charge_exception=?)
+                              bint raise_charge_exception=?, bint check_consistency=?)
 
     cpdef from_xyz(self, np.ndarray atomic_nums, np.ndarray coordinates, float critical_distance_factor=?, bint raise_atomtype_exception=?)
     
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index c628a4dd2e..d30d5bae8f 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -88,6 +88,8 @@ class Atom(Vertex):
     `mass`               ``int``             atomic mass of element (read only)
     `number`             ``int``             atomic number of element (read only)
     `symbol`             ``str``             atomic symbol of element (read only)
+    `site`               ``str``             type of adsorption site
+    `morphology`         ``str``             morphology of the adsorption site
     ==================== =================== ====================================
 
     Additionally, the ``mass``, ``number``, and ``symbol`` attributes of the
@@ -95,8 +97,8 @@ class Atom(Vertex):
     e.g. ``atom.symbol`` instead of ``atom.element.symbol``.
     """
 
-    def __init__(self, element=None, radical_electrons=0, charge=0, label='', lone_pairs=-100, coords=np.array([]),
-                 id=-1, props=None):
+    def __init__(self, element=None, radical_electrons=0, charge=0, label='', lone_pairs=-100, site='', morphology='', 
+                 coords=np.array([]), id=-1, props=None):
         Vertex.__init__(self)
         if isinstance(element, str):
             self.element = elements.__dict__[element]
@@ -107,6 +109,8 @@ def __init__(self, element=None, radical_electrons=0, charge=0, label='', lone_p
         self.label = label
         self.atomtype = None
         self.lone_pairs = lone_pairs
+        self.site = site 
+        self.morphology = morphology
         self.coords = coords
         self.id = id
         self.props = props or {}
@@ -139,6 +143,8 @@ def __reduce__(self):
             'sorting_label': self.sorting_label,
             'atomtype': self.atomtype.label if self.atomtype else None,
             'lone_pairs': self.lone_pairs,
+            'site': self.site,
+            'morphology': self.morphology,
         }
         if self.element.isotope == -1:
             element2pickle = self.element.symbol
@@ -157,6 +163,8 @@ def __setstate__(self, d):
         self.sorting_label = d['sorting_label']
         self.atomtype = ATOMTYPES[d['atomtype']] if d['atomtype'] else None
         self.lone_pairs = d['lone_pairs']
+        self.site = d['site']
+        self.morphology = d['morphology']
 
     def __hash__(self):
         """
@@ -250,6 +258,16 @@ def equivalent(self, other, strict=True):
                     if self.charge == charge: break
                 else:
                     return False
+            if ap.site:
+                for site in ap.site:
+                    if self.site == site: break
+                else:
+                    return False
+            if ap.morphology:
+                for morphology in ap.morphology:
+                    if self.morphology == morphology: break
+                else:
+                    return False
             if 'inRing' in self.props and 'inRing' in ap.props:
                 if self.props['inRing'] != ap.props['inRing']:
                     return False
@@ -293,6 +311,16 @@ def is_specific_case_of(self, other):
                         break
                 else:
                     return False
+            if atom.site:
+                for site in atom.site:
+                    if self.site == site: break
+                else:
+                    return False
+            if atom.morphology:
+                for morphology in atom.morphology:
+                    if self.morphology == morphology: break
+                else:
+                    return False
             if 'inRing' in self.props and 'inRing' in atom.props:
                 if self.props['inRing'] != atom.props['inRing']:
                     return False
@@ -316,6 +344,8 @@ def copy(self):
         a.label = self.label
         a.atomtype = self.atomtype
         a.lone_pairs = self.lone_pairs
+        a.site = self.site
+        a.morphology = self.morphology
         a.coords = self.coords[:]
         a.id = self.id
         a.props = deepcopy(self.props)
@@ -701,6 +731,8 @@ def get_order_str(self):
             return 'vdW'
         elif self.is_hydrogen_bond():
             return 'H'
+        elif self.is_reaction_bond():
+            return 'R'
         else:
             raise ValueError("Bond order {} does not have string representation.".format(self.order))
 
@@ -722,6 +754,8 @@ def set_order_str(self, new_order):
             self.order = 0
         elif new_order == 'H':
             self.order = 0.1
+        elif new_order == 'R':
+            self.order = 0.05
         else:
             # try to see if an float disguised as a string was input by mistake
             try:
@@ -815,6 +849,13 @@ def is_hydrogen_bond(self):
         """
         return self.is_order(0.1)
 
+    def is_reaction_bond(self):
+        """
+        Return ``True`` if the bond represents a reaction bond or ``False`` if
+        not.
+        """
+        return self.is_order(0.05)
+    
     def increment_order(self):
         """
         Update the bond as a result of applying a CHANGE_BOND action to
@@ -873,7 +914,7 @@ def get_bond_string(self):
         the atom labels in alphabetical order (i.e. 'C-H' is possible but not 'H-C')
         :return: str
         """
-        bond_symbol_mapping = {0.1: '~', 1: '-', 1.5: ':', 2: '=', 3: '#'}
+        bond_symbol_mapping = {0.05: '~', 0.1: '~', 1: '-', 1.5: ':', 2: '=', 3: '#'}
         atom_labels = [self.atom1.symbol, self.atom2.symbol]
         atom_labels.sort()
         try:
@@ -912,13 +953,16 @@ class Molecule(Graph):
     `inchi`                 ``str``     A string representation of the molecule in InChI
     `smiles`                ``str``     A string representation of the molecule in SMILES
     `fingerprint`           ``str``     A representation for fast comparison, set as molecular formula
+    `metal`                 ``str``     The metal of the metal surface the molecule is associated with
+    `facet`                 ``str``     The facet of the metal surface the molecule is associated with
     ======================= =========== ========================================
 
     A new molecule object can be easily instantiated by passing the `smiles` or
     `inchi` string representing the molecular structure.
     """
 
-    def __init__(self, atoms=None, symmetry=-1, multiplicity=-187, reactive=True, props=None, inchi='', smiles=''):
+    def __init__(self, atoms=None, symmetry=-1, multiplicity=-187, reactive=True, props=None, inchi='', smiles='', 
+                 metal='', facet=''):
         Graph.__init__(self, atoms)
         self.symmetry_number = symmetry
         self.multiplicity = multiplicity
@@ -927,6 +971,8 @@ def __init__(self, atoms=None, symmetry=-1, multiplicity=-187, reactive=True, pr
         self._inchi = None
         self._smiles = None
         self.props = props or {}
+        self.metal = metal
+        self.facet = facet
 
         if inchi and smiles:
             logging.warning('Both InChI and SMILES provided for Molecule instantiation, '
@@ -1001,7 +1047,7 @@ def __reduce__(self):
         """
         A helper function used when pickling an object.
         """
-        return (Molecule, (self.vertices, self.symmetry_number, self.multiplicity, self.reactive, self.props))
+        return (Molecule, (self.vertices, self.symmetry_number, self.multiplicity, self.reactive, self.props, self.metal, self.facet))
 
     @property
     def atoms(self):
@@ -1283,6 +1329,8 @@ def copy(self, deep=False):
             v2.sorting_label = v1.sorting_label
         other.multiplicity = self.multiplicity
         other.reactive = self.reactive
+        other.metal = self.metal 
+        other.facet = self.facet
         return other
 
     def merge(self, other):
@@ -1496,7 +1544,12 @@ def is_isomorphic(self, other, initial_map=None, generate_initial_map=False, sav
         # check multiplicity
         if self.multiplicity != other.multiplicity:
             return False
-        
+        #check metal
+        if self.metal != other.metal:
+            return False 
+        #check facet
+        if self.facet != other.facet:
+            return False
         # if given an initial map, ensure that it's valid.
         if initial_map:
             if not self.is_mapping_valid(other, initial_map, equivalent=True):
@@ -1534,7 +1587,12 @@ def find_isomorphism(self, other, initial_map=None, save_order=False, strict=Tru
         # check multiplicity
         if self.multiplicity != other.multiplicity:
             return []
-
+        #check metal
+        if self.metal != other.metal:
+            return []
+        #check facet
+        if self.facet != other.facet:
+            return []
         # Do the isomorphism comparison
         result = Graph.find_isomorphism(self, other, initial_map, save_order=save_order, strict=strict)
         return result
@@ -1561,7 +1619,12 @@ def is_subgraph_isomorphic(self, other, initial_map=None, generate_initial_map=F
         # Check multiplicity
         if group.multiplicity:
             if self.multiplicity not in group.multiplicity: return False
-
+        #check metal
+        if group.metal:
+            if self.metal not in group.metal: return False
+        #check facet
+        if group.facet:
+            if self.facet not in group.facet: return False
         # Compare radical counts
         if self.get_radical_count() < group.radicalCount:
             return False
@@ -1632,7 +1695,12 @@ def find_subgraph_isomorphisms(self, other, initial_map=None, save_order=False):
         # Check multiplicity
         if group.multiplicity:
             if self.multiplicity not in group.multiplicity: return []
-
+        #check metal
+        if group.metal:
+            if self.metal not in group.metal: return False
+        #check facet
+        if group.facet:
+            if self.facet not in group.facet: return False
         # Compare radical counts
         if self.get_radical_count() < group.radicalCount:
             return []
@@ -1718,7 +1786,7 @@ def from_smarts(self, smartsstr, raise_atomtype_exception=True):
         return self
 
     def from_adjacency_list(self, adjlist, saturate_h=False, raise_atomtype_exception=True,
-                            raise_charge_exception=True):
+                            raise_charge_exception=True, check_consistency=True):
         """
         Convert a string adjacency list `adjlist` to a molecular structure.
         Skips the first line (assuming it's a label) unless `withLabel` is
@@ -1726,7 +1794,8 @@ def from_adjacency_list(self, adjlist, saturate_h=False, raise_atomtype_exceptio
         """
         from rmgpy.molecule.adjlist import from_adjacency_list
 
-        self.vertices, self.multiplicity = from_adjacency_list(adjlist, group=False, saturate_h=saturate_h)
+        self.vertices, self.multiplicity, self.metal, self.facet = from_adjacency_list(adjlist, group=False, saturate_h=saturate_h,
+                                                               check_consistency=check_consistency)
         self.update_atomtypes(raise_exception=raise_atomtype_exception)
         self.identify_ring_membership()
 
@@ -1876,7 +1945,8 @@ def to_adjacency_list(self, label='', remove_h=False, remove_lone_pairs=False, o
         Convert the molecular structure to a string adjacency list.
         """
         from rmgpy.molecule.adjlist import to_adjacency_list
-        result = to_adjacency_list(self.vertices, self.multiplicity, label=label, group=False, remove_h=remove_h,
+        result = to_adjacency_list(self.vertices, self.multiplicity, metal=self.metal, facet=self.facet, 
+                                   label=label, group=False, remove_h=remove_h,
                                    remove_lone_pairs=remove_lone_pairs, old_style=old_style)
         return result
 
@@ -2296,7 +2366,8 @@ def to_group(self):
                                              label=atom.label,
                                              )
 
-        group = gr.Group(atoms=list(group_atoms.values()), multiplicity=[self.multiplicity])
+        group = gr.Group(atoms=list(group_atoms.values()), multiplicity=[self.multiplicity], metal=[self.metal],
+                         facet=[self.facet])
 
         # Create GroupBond for each bond between atoms in the molecule
         for atom in self.atoms: