Use rdkit for SSSR and RCs (bug fix + Python upgrade)

.conda/meta.yaml

@@ -64,7 +64,6 @@ requirements:
     - conda-forge::gprof2dot
     - conda-forge::numdifftools
     - conda-forge::quantities !=0.16.0,!=0.16.1
-    - conda-forge::ringdecomposerlib-python
     - rmg::pydas >=1.0.3
     - rmg::pydqed >=1.0.3
     - rmg::symmetry
@@ -114,7 +113,6 @@ requirements:
     - conda-forge::gprof2dot
     - conda-forge::numdifftools
     - conda-forge::quantities !=0.16.0,!=0.16.1
-    - conda-forge::ringdecomposerlib-python
     - rmg::pydas >=1.0.3
     - rmg::pydqed >=1.0.3
     - rmg::symmetry
@@ -165,7 +163,6 @@ test:
     - conda-forge::gprof2dot
     - conda-forge::numdifftools
     - conda-forge::quantities !=0.16.0,!=0.16.1
-    - conda-forge::ringdecomposerlib-python
     - rmg::pydas >=1.0.3
     - rmg::pydqed >=1.0.3
     - rmg::symmetry

.github/workflows/CI.yml

@@ -63,7 +63,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.9"]
+        python-version: ["3.9", "3.10", "3.11"]
         os: [macos-13, macos-latest, ubuntu-latest]
         include-rms: ["", "with RMS"]
         exclude:

.github/workflows/conda_build.yml

@@ -17,7 +17,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macos-13, macos-latest]
         numpy-version: ["1.26"]
-        python-version: ["3.9"]
+        python-version: ["3.9", "3.10", "3.11"]
     runs-on: ${{ matrix.os }}
     name: Build ${{ matrix.os }} Python ${{ matrix.python-version }} Numpy ${{ matrix.numpy-version }}
     defaults:

environment.yml

@@ -84,7 +84,6 @@ dependencies:
   # bug in quantities, see:
   # https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2694#issuecomment-2489286263
   - conda-forge::quantities !=0.16.0,!=0.16.1
-  - conda-forge::ringdecomposerlib-python
 
 # packages we maintain
   - rmg::pydas >=1.0.3

rmgpy/data/thermo.py

@@ -467,7 +467,7 @@ def is_ring_partial_matched(ring, matched_group):
     else:
         submol_ring, _ = convert_ring_to_sub_molecule(ring)
         sssr = submol_ring.get_smallest_set_of_smallest_rings()
-        sssr_grp = matched_group.get_smallest_set_of_smallest_rings()
+        sssr_grp = matched_group.make_sample_molecule().get_smallest_set_of_smallest_rings()
         if sorted([len(sr) for sr in sssr]) == sorted([len(sr_grp) for sr_grp in sssr_grp]):
             return False
         else:

rmgpy/molecule/adjlist.py

@@ -1093,7 +1093,7 @@ def to_adjacency_list(atoms, multiplicity, metal='', facet='', label=None, group
                 # numbers if doesn't work
                 try:
                     adjlist += bond.get_order_str()
-                except ValueError:
+                except (ValueError, TypeError):
                     adjlist += str(bond.get_order_num())
             adjlist += '}'
 

rmgpy/molecule/converter.pxd

@@ -29,7 +29,7 @@ cimport rmgpy.molecule.molecule as mm
 cimport rmgpy.molecule.element as elements
 
 
-cpdef to_rdkit_mol(mm.Molecule mol, bint remove_h=*, bint return_mapping=*, bint sanitize=*, bint save_order=?)
+cpdef to_rdkit_mol(mm.Molecule mol, bint remove_h=*, bint return_mapping=*, object sanitize=*, bint save_order=?)
 
 cpdef mm.Molecule from_rdkit_mol(mm.Molecule mol, object rdkitmol, bint raise_atomtype_exception=?)
 

rmgpy/molecule/converter.py

@@ -38,6 +38,7 @@
 import cython
 # Assume that rdkit is installed
 from rdkit import Chem
+from rdkit.Chem.rdchem import KekulizeException, AtomKekulizeException
 # Test if openbabel is installed
 try:
     from openbabel import openbabel
@@ -70,7 +71,9 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_o
     for index, atom in enumerate(mol.vertices):
         if atom.element.symbol == 'X':
             rd_atom = Chem.rdchem.Atom('Pt')  # not sure how to do this with linear scaling when this might not be Pt
-        else:
+        elif atom.element.symbol in ['R', 'L']:
+            rd_atom = Chem.rdchem.Atom(0)
+        else:   
             rd_atom = Chem.rdchem.Atom(atom.element.symbol)
         if atom.element.isotope != -1:
             rd_atom.SetIsotope(atom.element.isotope)
@@ -96,8 +99,9 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_o
             else:
                 label_dict[label] = [saved_index]
     rd_bonds = Chem.rdchem.BondType
+    # no vdW bond in RDKit, so "ZERO" or "OTHER" might be OK
     orders = {'S': rd_bonds.SINGLE, 'D': rd_bonds.DOUBLE, 'T': rd_bonds.TRIPLE, 'B': rd_bonds.AROMATIC,
-              'Q': rd_bonds.QUADRUPLE}
+              'Q': rd_bonds.QUADRUPLE, 'vdW': rd_bonds.ZERO, 'H': rd_bonds.HYDROGEN, 'R': rd_bonds.UNSPECIFIED} 
     # Add the bonds
     for atom1 in mol.vertices:
         for atom2, bond in atom1.edges.items():
@@ -107,7 +111,10 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_o
             index2 = atoms.index(atom2)
             if index1 < index2:
                 order_string = bond.get_order_str()
-                order = orders[order_string]
+                if order_string is None:
+                    order = rd_bonds.UNSPECIFIED
+                else:
+                    order = orders[order_string]
                 rdkitmol.AddBond(index1, index2, order)
 
     # Make editable mol into a mol and rectify the molecule
@@ -119,8 +126,14 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_o
     for atom in rdkitmol.GetAtoms():
         if atom.GetAtomicNum() > 1:
             atom.SetNoImplicit(True)
-    if sanitize:
+    if sanitize == True:
         Chem.SanitizeMol(rdkitmol)
+    elif sanitize == "partial":
+        try:
+            Chem.SanitizeMol(rdkitmol, sanitizeOps=Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES)
+        except (KekulizeException, AtomKekulizeException):
+            logging.debug("Kekulization failed; sanitizing without Kekulize")
+            Chem.SanitizeMol(rdkitmol, sanitizeOps=Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES ^ Chem.SANITIZE_KEKULIZE)
     if remove_h:
         rdkitmol = Chem.RemoveHs(rdkitmol, sanitize=sanitize)
     if return_mapping:

rmgpy/molecule/draw.py

@@ -155,7 +155,7 @@ def clear(self):
         self.surface = None
         self.cr = None
 
-    def draw(self, molecule, file_format, target=None):
+    def draw(self, molecule, file_format, target=None, use_rdkit=True):
         """
         Draw the given `molecule` using the given image `file_format` - pdf, svg, ps, or
         png. If `path` is given, the drawing is saved to that location on disk. The
@@ -165,6 +165,9 @@ def draw(self, molecule, file_format, target=None):
         This function returns the Cairo surface and context used to create the
         drawing, as well as a bounding box for the molecule being drawn as the
         tuple (`left`, `top`, `width`, `height`).
+
+        If `use_rdkit` is True, then the RDKit 2D coordinate generation is used to generate the coordinates.
+        If `use_rdkit` is False, then the molecule is drawn using our (deprecated) original algorithm.
         """
 
         # The Cairo 2D graphics library (and its Python wrapper) is required for
@@ -219,13 +222,13 @@ def draw(self, molecule, file_format, target=None):
                 if molecule.contains_surface_site():
                     try:
                         self._connect_surface_sites()
-                        self._generate_coordinates()
+                        self._generate_coordinates(use_rdkit=use_rdkit)
                         self._disconnect_surface_sites()
                     except AdsorbateDrawingError as e:
                         self._disconnect_surface_sites()
-                        self._generate_coordinates(fix_surface_sites=False)
+                        self._generate_coordinates(fix_surface_sites=False, use_rdkit=use_rdkit)
                 else:
-                    self._generate_coordinates()
+                    self._generate_coordinates(use_rdkit=use_rdkit)
                 self._replace_bonds(old_bond_dictionary)
 
                 # Generate labels to use
@@ -341,7 +344,7 @@ def _find_ring_groups(self):
                 if not found:
                     self.ringSystems.append([cycle])
 
-    def _generate_coordinates(self, fix_surface_sites=True):
+    def _generate_coordinates(self, fix_surface_sites=True, use_rdkit=True):
         """
         Generate the 2D coordinates to be used when drawing the current
         molecule. The function uses rdKits 2D coordinate generation.
@@ -372,15 +375,34 @@ def _generate_coordinates(self, fix_surface_sites=True):
                 self.coordinates[1, :] = [0.5, 0.0]
             return self.coordinates
 
-        # Decide whether we can use RDKit or have to generate coordinates ourselves
-        for atom in self.molecule.atoms:
-            if atom.charge != 0:
-                use_rdkit = False
-                break
-        else:  # didn't break
-            use_rdkit = True
+        if use_rdkit == True:
+            # Use RDKit 2D coordinate generation:
+
+            # Generate the RDkit molecule from the RDkit molecule, use geometry
+            # in order to match the atoms in the rdmol with the atoms in the
+            # RMG molecule (which is required to extract coordinates).
+            self.geometry = Geometry(None, None, self.molecule, None)
+
+            rdmol, rd_atom_idx = self.geometry.rd_build()
+            AllChem.Compute2DCoords(rdmol)
+
+            # Extract the coordinates from each atom.
+            for atom in atoms:
+                index = rd_atom_idx[atom]
+                point = rdmol.GetConformer(0).GetAtomPosition(index)
+                coordinates[index, :] = [point.x * 0.6, point.y * 0.6]
+
+            # RDKit generates some molecules more vertically than horizontally,
+            # Especially linear ones. This will reflect any molecule taller than
+            # it is wide across the line y=x
+            ranges = np.ptp(coordinates, axis=0)
+            if ranges[1] > ranges[0]:
+                temp = np.copy(coordinates)
+                coordinates[:, 0] = temp[:, 1]
+                coordinates[:, 1] = temp[:, 0]
 
-        if not use_rdkit:
+        else:
+            logging.warning("Using deprecated molecule drawing algorithm; undesired behavior may occur. Consider using use_rdkit=True.")
             if len(self.cycles) > 0:
                 # Cyclic molecule
                 backbone = self._find_cyclic_backbone()
@@ -438,32 +460,6 @@ def _generate_coordinates(self, fix_surface_sites=True):
             # minimize likelihood of overlap
             self._generate_neighbor_coordinates(backbone)
 
-        else:
-            # Use RDKit 2D coordinate generation:
-
-            # Generate the RDkit molecule from the RDkit molecule, use geometry
-            # in order to match the atoms in the rdmol with the atoms in the
-            # RMG molecule (which is required to extract coordinates).
-            self.geometry = Geometry(None, None, self.molecule, None)
-
-            rdmol, rd_atom_idx = self.geometry.rd_build()
-            AllChem.Compute2DCoords(rdmol)
-
-            # Extract the coordinates from each atom.
-            for atom in atoms:
-                index = rd_atom_idx[atom]
-                point = rdmol.GetConformer(0).GetAtomPosition(index)
-                coordinates[index, :] = [point.x * 0.6, point.y * 0.6]
-
-            # RDKit generates some molecules more vertically than horizontally,
-            # Especially linear ones. This will reflect any molecule taller than
-            # it is wide across the line y=x
-            ranges = np.ptp(coordinates, axis=0)
-            if ranges[1] > ranges[0]:
-                temp = np.copy(coordinates)
-                coordinates[:, 0] = temp[:, 1]
-                coordinates[:, 1] = temp[:, 0]
-
         # For surface species
         if fix_surface_sites and self.molecule.contains_surface_site():
             if len(self.molecule.atoms) == 1:

rmgpy/molecule/fragment.py

@@ -490,24 +490,30 @@ def get_formula(self):
 
         return formula
 
-    def to_rdkit_mol(self, remove_h=False, return_mapping=True, save_order=False):
+    def to_rdkit_mol(self, *args, **kwargs):
         """
         Convert a molecular structure to a RDKit rdmol object.
         """
-        if remove_h:
+        remove_h = kwargs.get('remove_h', False)
+
+        if remove_h == True:
             # because we're replacing
             # cutting labels with hydrogens
             # so do not allow removeHs to be True
             raise "Currently fragment to_rdkit_mol only allows keeping all the hydrogens."
 
         mol0, mapping = self.get_representative_molecule("minimal", update=False)
 
+        return_mapping = kwargs.get("return_mapping", False)
+        if return_mapping == False:
+            kwargs["return_mapping"] = True
+            logging.warning("Fragment to_rdkit_mol expects to return a tuple. "
+                            "Setting return_mapping = True; please double-check your code to ensure this is what you want.")
+
         rdmol, rdAtomIdx_mol0 = converter.to_rdkit_mol(
             mol0,
-            remove_h=remove_h,
-            return_mapping=return_mapping,
-            sanitize=True,
-            save_order=save_order,
+            *args,
+            **kwargs
         )
 
         rdAtomIdx_frag = {}
@@ -587,7 +593,7 @@ def get_aromatic_rings(self, rings=None, save_order=False):
         """
         Returns all aromatic rings as a list of atoms and a list of bonds.
 
-        Identifies rings using `Graph.get_smallest_set_of_smallest_rings()`, then uses RDKit to perceive aromaticity.
+        Identifies rings, then uses RDKit to perceive aromaticity.
         RDKit uses an atom-based pi-electron counting algorithm to check aromaticity based on Huckel's Rule.
         Therefore, this method identifies "true" aromaticity, rather than simply the RMG bond type.
 
@@ -911,7 +917,7 @@ def sliceitup_arom(self, molecule, size_threshold=None):
         _, cutting_label_list = self.detect_cutting_label(molecule_smiles)
         # transfer to rdkit molecule for substruct matching
         f = self.from_smiles_like_string(molecule_smiles)
-        molecule_to_cut, rdAtomIdx_frag = f.to_rdkit_mol()
+        molecule_to_cut, rdAtomIdx_frag = f.to_rdkit_mol(remove_h=False, return_mapping=True)
 
         # replace CuttingLabel to special Atom (metal) in rdkit
         for atom, idx in rdAtomIdx_frag.items():
@@ -1037,7 +1043,7 @@ def sliceitup_aliph(self, molecule, size_threshold=None):
         _, cutting_label_list = self.detect_cutting_label(molecule_smiles)
         # transfer to rdkit molecule for substruct matching
         f = self.from_smiles_like_string(molecule_smiles)
-        molecule_to_cut, rdAtomIdx_frag = f.to_rdkit_mol()
+        molecule_to_cut, rdAtomIdx_frag = f.to_rdkit_mol(remove_h=False, return_mapping=True)
 
         # replace CuttingLabel to special Atom (metal) in rdkit
         for atom, idx in rdAtomIdx_frag.items():

rmgpy/molecule/graph.pxd

@@ -127,16 +127,6 @@ cdef class Graph(object):
     cpdef bint _is_chain_in_cycle(self, list chain) except -2
 
     cpdef list get_all_cyclic_vertices(self)
-
-    cpdef list get_all_polycyclic_vertices(self)
-
-    cpdef list get_polycycles(self)
-
-    cpdef list get_monocycles(self)
-
-    cpdef tuple get_disparate_cycles(self)
-
-    cpdef tuple _merge_cycles(self, list cycle_sets)
 
     cpdef list get_all_cycles(self, Vertex starting_vertex)
 
@@ -146,13 +136,7 @@ cdef class Graph(object):
 
     cpdef list _explore_cycles_recursively(self, list chain, list cycles)
 
-    cpdef list get_smallest_set_of_smallest_rings(self)
-
-    cpdef list get_relevant_cycles(self)
-
     cpdef list sort_cyclic_vertices(self, list vertices)
-
-    cpdef int get_max_cycle_overlap(self)
 
     cpdef list get_largest_ring(self, Vertex vertex)