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,9 @@ 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):
+            adjlist += ' n{0}'.format(atom_ncoord[atom])
         # Isotopes
         if atom_isotope[atom] != -1:
             adjlist += ' i{0}'.format(atom_isotope[atom])

rmgpy/molecule/converter.py

@@ -161,6 +161,7 @@ def from_rdkit_mol(mol, rdkitmol, raise_atomtype_exception=True):
                    bond=mm.Bond)
 
     mol.vertices = []
+    mol.invalidate_cycle_cache()
 
     # Add hydrogen atoms to complete molecule if needed
     rdkitmol.UpdatePropertyCache(strict=False)
@@ -314,6 +315,7 @@ def from_ob_mol(mol, obmol, raise_atomtype_exception=True):
         raise DependencyError('OpenBabel is not installed. Please install or use RDKit.')
 
     mol.vertices = []
+    mol.invalidate_cycle_cache()
 
     # Add hydrogen atoms to complete molecule if needed
     obmol.AddHydrogens()

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):
@@ -672,6 +675,7 @@ def to_smiles(self):
                 final_vertices.append(atom)
 
         smiles_before.vertices = final_vertices
+        smiles_before.invalidate_cycle_cache()
         mol_repr = Molecule()
         mol_repr.atoms = smiles_before.vertices
         mol_repr.update()
@@ -786,6 +790,7 @@ def from_rdkit_mol(self, rdkitmol, atom_replace_dict=None):
         from rdkit import Chem
 
         self.vertices = []
+        self.invalidate_cycle_cache()
 
         # Add hydrogen atoms to complete molecule if needed
         rdkitmol.UpdatePropertyCache(strict=False)

rmgpy/molecule/graph.pxd

@@ -37,12 +37,13 @@ cdef class Vertex(object):
     cdef public bint terminal
     cdef public Vertex mapping
     cdef public bint ignore
-
+    cdef public bint excluded
+
     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)
 
@@ -73,9 +74,11 @@ cdef Vertex _get_edge_vertex2(Edge edge)
 cdef class Graph(object):
 
     cdef public list vertices
-    
+
     cdef public list ordered_vertices
 
+    cdef list _relevant_cycles
+
     cpdef Vertex add_vertex(self, Vertex vertex)
 
     cpdef Edge add_edge(self, Edge edge)
@@ -110,36 +113,46 @@ 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=?, bint check_labels=?)
+
+    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?) except -2
 
-    cpdef list find_isomorphism(self, Graph other, dict initial_map=?, bint save_order=?, bint strict=?)
+    cpdef list find_subgraph_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?)
 
-    cpdef bint is_subgraph_isomorphic(self, Graph other, dict initial_map=?, bint save_order=?) except -2
+    cpdef bint is_intersection_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_intersection_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?)
+
+    cpdef list find_largest_incomplete_isomorphisms(self, Graph other, dict initial_map=?, bint save_order=?, bint check_labels=?, bint find_all=?)
 
     cpdef bint is_cyclic(self) except -2
 
     cpdef bint is_vertex_in_cycle(self, Vertex vertex) except -2
 
     cpdef bint is_edge_in_cycle(self, Edge edge) except -2
 
-    cpdef bint _is_chain_in_cycle(self, list chain) except -2
+    cdef list _get_relevant_cycles(self)
+
+    cpdef invalidate_cycle_cache(self)
 
     cpdef list get_all_cyclic_vertices(self)
 
     cpdef list get_all_cycles(self, Vertex starting_vertex)
 
+    cpdef list _explore_cycles_recursively(self, list chain, list cycles)
+
     cpdef list get_all_cycles_of_size(self, int size)
 
     cpdef list get_all_simple_cycles_of_size(self, int size)
 
-    cpdef list _explore_cycles_recursively(self, list chain, list cycles)
-
     cpdef list sort_cyclic_vertices(self, list vertices)
 
     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=?, bint intersection=?) 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=?)