Isomorphism Improvements

.github/workflows/CI.yml

@@ -37,7 +37,7 @@ concurrency:
 env:
   # if running on RMG-Py but requiring changes on an un-merged branch of RMG-database, replace
   # main with the name of the branch
-  RMG_DATABASE_BRANCH: main
+  RMG_DATABASE_BRANCH: atomtype_renaming
   # RMS branch to use for ReactionMechanismSimulator installation
   RMS_BRANCH: for_rmg
   # RMS mode used for install_rms.sh

documentation/source/reference/molecule/atomtype.rst

@@ -30,10 +30,10 @@ Atom type       Description
 *Carbon atom types*
 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 ``C``           carbon atom with any local bond structure
-``Ca``          carbon atom with two lone pairs and no bonds
-``Cs``          carbon atom with up to four single bonds
+``Car``         carbon atom with two lone pairs and no bonds
+``Css``         carbon atom with up to four single bonds
 ``Csc``         charged carbon atom with up to three single bonds
-``Cd``          carbon atom with one double bond (not to O or S) and up to two single bonds
+``Cdb``         carbon atom with one double bond (not to O or S) and up to two single bonds
 ``Cdc``         charged carbon atom with one double bond and up to one single bond
 ``CO``          carbon atom with one double bond to oxygen and up to two single bonds
 ``CS``          carbon atom with one double bond to sulfur and up to two single bonds

rmgpy/molecule/adjlist.py

@@ -753,7 +753,21 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False, check_consistenc
                 for m in m_state:
                     morphologies.append(m[1:-1])
                 index += 1
-
+
+        props = {}
+
+        # Next the coordination numbers (if provided)
+        if len(data) > index:
+            n_state = data[index]
+            if n_state[0] == 'n':
+                if n_state[1] == '[':
+                    n_state = n_state[2:-1].split(',')
+                else:
+                    n_state = [n_state[1:]]
+                coords = [int(n) for n in n_state]
+                props['Ncoord'] = coords
+                index += 1
+
         # Next the isotope (if provided)
         isotope = -1
         if len(data) > index:
@@ -763,7 +777,6 @@ def from_adjacency_list(adjlist, group=False, saturate_h=False, check_consistenc
                 index += 1
 
         # Next ring membership info (if provided)
-        props = {}
         if len(data) > index:
             r_state = data[index]
             if r_state[0] == 'r':
@@ -943,6 +956,7 @@ def to_adjacency_list(atoms, multiplicity, metal='', facet='', label=None, group
     atom_props = {}
     atom_site = {}
     atom_morphology = {}
+    atom_ncoord = {}
     if group:
         for atom in atom_numbers:
             # Atom type(s)
@@ -989,7 +1003,17 @@ def to_adjacency_list(atoms, multiplicity, metal='', facet='', label=None, group
                 atom_morphology[atom] = None  # Empty list indicates wildcard
             else:
                 atom_morphology[atom] = '["{0}"]'.format('","'.join(s for s in atom.morphology))
-
+
+            # Coordination Number
+            if 'Ncoord' not in atom.props:
+                atom_ncoord[atom] = None
+            elif len(atom.props['Ncoord']) == 1:
+                atom_ncoord[atom] = atom.props['Ncoord'][0]
+            elif len(atom.props['Ncoord']) == 0:
+                atom_ncoord[atom] = []
+            else:
+                atom_ncoord[atom] = '[{0}]'.format(','.join(str(s) for s in atom.props['Ncoord']))
+
             # Isotopes
             atom_isotope[atom] = -1
 
@@ -1058,6 +1082,10 @@ def to_adjacency_list(atoms, multiplicity, metal='', facet='', label=None, group
         # Morphologies
         if atom_morphology[atom]:
             adjlist += ' m{0}'.format(atom_morphology[atom])
+        # Coordination numbers
+        if group and atom_ncoord[atom] is not None and (isinstance(atom_ncoord[atom], int) or len(atom_ncoord[atom]) > 0):
+            if atom_ncoord[atom]:
+                adjlist += ' n{0}'.format(atom_ncoord[atom])
         # Isotopes
         if atom_isotope[atom] != -1:
             adjlist += ' i{0}'.format(atom_isotope[atom])

rmgpy/molecule/fragment.py

@@ -80,20 +80,23 @@ def __repr__(self):
     def symbol(self):
         return self.name
 
-    def is_specific_case_of(self, other):
+    def is_specific_case_of(self, other, check_labels=False):
         """
         Return ``True`` if `self` is a specific case of `other`, or ``False``
         otherwise. At this moment, this is the same as the :math:`equivalent()`.
         """
-        return self.equivalent(other)
+        return self.equivalent(other, check_labels=check_labels)
 
-    def equivalent(self, other, strict=True):
+    def equivalent(self, other, strict=True, check_labels=False):
         """
         Return ``True`` if `other` is indistinguishable from this CuttingLabel, or
         ``False`` otherwise. If `other` is an :class:`CuttingLabel` object, then all
-        attributes must match exactly.
+        attributes must match exactly. If ``check_labels`` is ``True``, the
+        `label` attributes must also match.
         """
         if isinstance(other, CuttingLabel):
+            if check_labels and self.label != other.label:
+                return False
             return self.name == other.name
         else:
             return False
@@ -247,7 +250,7 @@ def get_radical_count(self):
         return radicals
 
     def is_subgraph_isomorphic(
-        self, other, initial_map=None, generate_initial_map=False, save_order=False
+        self, other, initial_map=None, generate_initial_map=False, save_order=False, check_labels=False
     ):
         """
         Fragment's subgraph isomorphism check is done by first creating
@@ -303,15 +306,15 @@ def is_subgraph_isomorphic(
                     for i, key in enumerate(keys):
                         initial_map[key] = atmlist[i]
                     if self.is_mapping_valid(
-                        other, initial_map, equivalent=False
+                        other, initial_map, equivalent=False, strict=True, check_labels=check_labels
                     ) and Graph.is_subgraph_isomorphic(
-                        self, other, initial_map, save_order=save_order
+                        self, other, initial_map, save_order=save_order, check_labels=check_labels
                     ):
                         return True
                 else:
                     return False
             else:
-                if not self.is_mapping_valid(other, initial_map, equivalent=False):
+                if not self.is_mapping_valid(other, initial_map, equivalent=False, strict=True, check_labels=check_labels):
                     return False
 
         # Do the isomorphism comparison
@@ -322,7 +325,7 @@ def is_subgraph_isomorphic(
                 repr_mol_vertex = mapping[fragment_vertex]
                 new_initial_map[repr_mol_vertex] = initial_map[fragment_vertex]
 
-        result = Graph.is_subgraph_isomorphic(self.mol_repr, other, new_initial_map)
+        result = Graph.is_subgraph_isomorphic(self.mol_repr, other, new_initial_map, check_labels=check_labels)
         return result
 
     def calculate_cp0(self):

rmgpy/molecule/graph.pxd

@@ -40,9 +40,9 @@ cdef class Vertex(object):
 
     cpdef Vertex copy(self)
 
-    cpdef bint equivalent(self, Vertex other, bint strict=?) except -2
+    cpdef bint equivalent(self, Vertex other, bint strict=?, bint check_labels=?) except -2
 
-    cpdef bint is_specific_case_of(self, Vertex other) except -2
+    cpdef bint is_specific_case_of(self, Vertex other, bint check_labels=?) except -2
 
     cpdef reset_connectivity_values(self)
 
@@ -110,13 +110,17 @@ cdef class Graph(object):
 
     cpdef restore_vertex_order(self)
 
-    cpdef bint is_isomorphic(self, Graph other, dict initial_map=?, bint generate_initial_map=?, bint save_order=?, bint strict=?) except -2
+    cpdef bint is_isomorphic(self, Graph other, dict initial_map=?, bint generate_initial_map=?, bint save_order=?, bint strict=?, bint check_labels=?) except -2
 
-    cpdef list find_isomorphism(self, Graph other, dict initial_map=?, bint save_order=?, bint strict=?)
+    cpdef list find_isomorphism(self, Graph other, dict initial_map=?, bint save_order=?, bint strict=?, bint check_labels=?)
 
-    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=?, bint save_order=?) except -2
+    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?) except -2
 
-    cpdef list find_subgraph_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?)
+    cpdef list find_subgraph_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?)
+
+    cpdef bint is_intersection_isomorphic(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?) except -2
+
+    cpdef list find_intersection_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?)
 
     cpdef bint is_cyclic(self) except -2
 
@@ -140,6 +144,8 @@ cdef class Graph(object):
 
     cpdef list get_largest_ring(self, Vertex vertex)
 
-    cpdef bint is_mapping_valid(self, Graph other, dict mapping, bint equivalent=?, bint strict=?) except -2
+    cpdef bint is_mapping_valid(self, Graph other, dict mapping, bint equivalent=?, bint strict=?, bint check_labels=?) except -2
+
+    cpdef bint has_same_labels(self, Graph other, list ignore_labels=?) except -2
 
     cpdef list get_edges_in_cycle(self, list vertices, bint sort=?)

rmgpy/molecule/graph.pyx

@@ -95,19 +95,23 @@ cdef class Vertex(object):
         new = Vertex()
         return new
 
-    cpdef bint equivalent(self, Vertex other, bint strict=True) except -2:
+    cpdef bint equivalent(self, Vertex other, bint strict=True, bint check_labels=False) except -2:
         """
         Return :data:`True` if two vertices `self` and `other` are semantically
         equivalent, or :data:`False` if not. You should reimplement this
         function in a derived class if your vertices have semantic information.
+        If `check_labels` is ``True``, subclasses with a `label` attribute
+        should also require that the labels match.
         """
         return True
 
-    cpdef bint is_specific_case_of(self, Vertex other) except -2:
+    cpdef bint is_specific_case_of(self, Vertex other, bint check_labels=False) except -2:
         """
         Return ``True`` if `self` is semantically more specific than `other`,
         or ``False`` if not. You should reimplement this function in a derived
-        class if your edges have semantic information.
+        class if your edges have semantic information. If `check_labels` is
+        ``True``, subclasses with a `label` attribute should also require that
+        the labels match.
         """
         return True
 
@@ -499,7 +503,7 @@ cdef class Graph(object):
         else:
             self.vertices = self.ordered_vertices
 
-    cpdef bint is_isomorphic(self, Graph other, dict initial_map=None, bint generate_initial_map=False, bint save_order=False, bint strict=True) except -2:
+    cpdef bint is_isomorphic(self, Graph other, dict initial_map=None, bint generate_initial_map=False, bint save_order=False, bint strict=True, bint check_labels=False) except -2:
         """
         Returns :data:`True` if two graphs are isomorphic and :data:`False`
         otherwise. Uses the VF2 algorithm of Vento and Foggia.
@@ -509,6 +513,7 @@ cdef class Graph(object):
             generate_initial_map (bool, optional): if ``True``, initialize map by pairing atoms with same labels
             save_order (bool, optional):  if ``True``, reset atom order after performing atom isomorphism
             strict (bool, optional):     if ``False``, perform isomorphism ignoring electrons
+            check_labels (bool, optional): if ``True``, atoms only match if their `label` attributes match
         """
         if generate_initial_map:
             initial_map = dict()
@@ -520,12 +525,12 @@ cdef class Graph(object):
                             break
                     else:
                         return False
-            if not self.is_mapping_valid(other, initial_map, equivalent=True):
+            if not self.is_mapping_valid(other, initial_map, equivalent=True, strict=True, check_labels=check_labels):
                 return False
 
-        return vf2.is_isomorphic(self, other, initial_map, save_order=save_order, strict=strict)
+        return vf2.is_isomorphic(self, other, initial_map, save_order=save_order, strict=strict, check_labels=check_labels)
 
-    cpdef list find_isomorphism(self, Graph other, dict initial_map=None, bint save_order=False, bint strict=True):
+    cpdef list find_isomorphism(self, Graph other, dict initial_map=None, bint save_order=False, bint strict=True, bint check_labels=False):
         """
         Returns :data:`True` if `other` is subgraph isomorphic and :data:`False`
         otherwise, and the matching mapping.
@@ -535,24 +540,41 @@ cdef class Graph(object):
             initial_map (dict, optional): initial atom mapping to use
             save_order (bool, optional):  if ``True``, reset atom order after performing atom isomorphism
             strict (bool, optional):     if ``False``, perform isomorphism ignoring electrons
+            check_labels (bool, optional): if ``True``, atoms only match if their `label` attributes match
         """
-        return vf2.find_isomorphism(self, other, initial_map, save_order=save_order, strict=strict)
+        return vf2.find_isomorphism(self, other, initial_map, save_order=save_order, strict=strict, check_labels=check_labels)
 
-    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=None, bint save_order=False) except -2:
+    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=None, bint save_order=False, bint check_labels=False) except -2:
         """
         Returns :data:`True` if `other` is subgraph isomorphic and :data:`False`
         otherwise. Uses the VF2 algorithm of Vento and Foggia.
         """
-        return vf2.is_subgraph_isomorphic(self, other, initial_map, save_order=save_order)
+        return vf2.is_subgraph_isomorphic(self, other, initial_map, save_order=save_order, check_labels=check_labels)
 
-    cpdef list find_subgraph_isomorphisms(self, Graph other, dict initial_map=None, bint save_order=False):
+    cpdef list find_subgraph_isomorphisms(self, Graph other, dict initial_map=None, bint save_order=False, bint check_labels=False):
         """
         Returns :data:`True` if `other` is subgraph isomorphic and :data:`False`
         otherwise. Also returns the lists all of valid mappings.
 
         Uses the VF2 algorithm of Vento and Foggia.
         """
-        return vf2.find_subgraph_isomorphisms(self, other, initial_map, save_order=save_order)
+        return vf2.find_subgraph_isomorphisms(self, other, initial_map, save_order=save_order, check_labels=check_labels)
+
+    cpdef bint is_intersection_isomorphic(self, Graph other, dict initial_map=None, bint save_order=False, bint check_labels=False) except -2:
+        """
+        Returns :data:`True` if `other` is intersection isomorphic and :data:`False`
+        otherwise. Uses the VF2 algorithm of Vento and Foggia.
+        """
+        return vf2.is_intersection_isomorphic(self, other, initial_map, save_order=save_order, check_labels=check_labels)
+
+    cpdef list find_intersection_isomorphisms(self, Graph other, dict initial_map=None, bint save_order=False, bint check_labels=False):
+        """
+        Returns :data:`True` if `other` is intersection isomorphic and :data:`False`
+        otherwise. Also returns the lists all of valid mappings.
+
+        Uses the VF2 algorithm of Vento and Foggia.
+        """
+        return vf2.find_intersection_isomorphisms(self, other, initial_map, save_order=save_order, check_labels=check_labels)
 
     cpdef bint is_cyclic(self) except -2:
         """
@@ -838,14 +860,16 @@ cdef class Graph(object):
                 longest_cycle = cycle
         return longest_cycle
 
-    cpdef bint is_mapping_valid(self, Graph other, dict mapping, bint equivalent=True, bint strict=True) except -2:
+    cpdef bint is_mapping_valid(self, Graph other, dict mapping, bint equivalent=True, bint strict=True, bint check_labels=False) except -2:
         """
         Check that a proposed `mapping` of vertices from `self` to `other`
         is valid by checking that the vertices and edges involved in the
         mapping are mutually equivalent.  If equivalent is ``True`` it checks
         if atoms and edges are equivalent, if ``False`` it checks if they
         are specific cases of each other. If strict is ``True``, electrons
-        and bond orders are considered, and ignored if ``False``.
+        and bond orders are considered, and ignored if ``False``. If
+        check_labels is ``True``, atoms only match if their `label`
+        attributes also match.
         """
         cdef Vertex vertex1, vertex2
         cdef list vertices1, vertices2
@@ -855,10 +879,10 @@ cdef class Graph(object):
         # Check that the mapped pairs of vertices compare True
         for vertex1, vertex2 in mapping.items():
             if equivalent:
-                if not vertex1.equivalent(vertex2, strict=strict):
+                if not vertex1.equivalent(vertex2, strict=strict, check_labels=check_labels):
                     return False
             else:
-                if not vertex1.is_specific_case_of(vertex2):
+                if not vertex1.is_specific_case_of(vertex2, check_labels=check_labels):
                     return False
 
         # Check that any edges connected mapped vertices are equivalent
@@ -890,6 +914,36 @@ cdef class Graph(object):
         # mapping is valid
         return True
 
+    cpdef bint has_same_labels(self, Graph other, list ignore_labels=None) except -2:
+        """
+        Returns ``True`` if `self` and `other` have the same labels on
+        their vertices, with the same number of vertices bearing each
+        label (i.e. the multisets of vertex labels are equal). Vertices
+        without a label (an empty or unset `label` attribute) are ignored.
+        Returns ``False`` otherwise.
+
+        If `ignore_labels` is given, vertices bearing any of those labels
+        are excluded from the comparison entirely.
+        """
+        cdef dict labels1, labels2
+        cdef set skip
+
+        skip = set(ignore_labels) if ignore_labels else set()
+
+        labels1 = {}
+        for vertex in self.vertices:
+            label = vertex.label
+            if label and label not in skip:
+                labels1[label] = labels1.get(label, 0) + 1
+
+        labels2 = {}
+        for vertex in other.vertices:
+            label = vertex.label
+            if label and label not in skip:
+                labels2[label] = labels2.get(label, 0) + 1
+
+        return labels1 == labels2
+
     cpdef list get_edges_in_cycle(self, list vertices, bint sort=False):
         """
         For a given list of atoms comprising a ring, return the set of bonds