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,12 +63,9 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.9"]
-        os: [macos-13, macos-latest, ubuntu-latest]
+        python-version: ["3.9", "3.10", "3.11"]
+        os: [macos-latest, ubuntu-latest]
         include-rms: ["", "with RMS"]
-        exclude:
-          - os: macos-13  # GitHub's runners just aren't up to the task of installing Julia
-            include-rms: 'with RMS'
     runs-on: ${{ matrix.os }}
     name: Python ${{ matrix.python-version }} ${{ matrix.os }} Build and Test ${{ matrix.include-rms }}
     # skip scheduled runs from forks

.github/workflows/conda_build.yml

@@ -15,9 +15,9 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, macos-13, macos-latest]
+        os: [ubuntu-latest, 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:

Makefile

@@ -28,7 +28,7 @@ clean-solver:
 
 install:
 	@ python utilities.py check-pydas
-	python -m pip install -vv -e .
+	python -m pip install --no-build-isolation -vv -e .
 
 q2dtor:
 	@ echo -e "\nInstalling Q2DTor...\n"

environment.yml

@@ -58,7 +58,7 @@ dependencies:
   - conda-forge::rdkit >=2022.09.1
 
 # Python tools
-  - conda-forge::python >=3.9  # leave as GEQ so that GitHub actions can add EQ w/o breaking (contradictory deps)
+  - conda-forge::python >=3.9,<3.12  # leave as GEQ so that GitHub actions can add EQ w/o breaking (contradictory deps)
   - conda-forge::setuptools <80
   - conda-forge::coverage
   - conda-forge::cython >=0.25.2,<3.1
@@ -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:
@@ -483,7 +483,7 @@ def bicyclic_decomposition_for_polyring(polyring):
     """
 
     submol, _ = convert_ring_to_sub_molecule(polyring)
-    sssr = submol.get_deterministic_sssr()
+    sssr = submol.get_smallest_set_of_smallest_rings()
 
     ring_pair_with_common_atoms_list = []
     ring_occurances_dict = {}
@@ -555,7 +555,7 @@ def split_bicyclic_into_single_rings(bicyclic_submol):
     Splits a given bicyclic submolecule into two individual single 
     ring submolecules (a list of `Molecule`s ).
     """
-    sssr = bicyclic_submol.get_deterministic_sssr()
+    sssr = bicyclic_submol.get_smallest_set_of_smallest_rings()
 
     return [convert_ring_to_sub_molecule(sssr[0])[0],
             convert_ring_to_sub_molecule(sssr[1])[0]]

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=?, bint ignore_bond_orders=?)
 
 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
@@ -49,55 +50,61 @@
 from rmgpy.exceptions import DependencyError
 
 
-def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_order=False):
+def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True,
+                 save_order=False, ignore_bond_orders=False):
     """
     Convert a molecular structure to a RDKit rdmol object. Uses
     `RDKit <https://rdkit.org/>`_ to perform the conversion.
     Perceives aromaticity and, unless remove_h==False, removes Hydrogen atoms.
 
     If return_mapping==True then it also returns a dictionary mapping the
     atoms to RDKit's atom indices.
+
+    If ignore_bond_orders==True, all bonds are converted to unknown bonds, and 
+    sanitization is skipped. This is helpful when all you want is ring perception,
+    for example. Must also set sanitize=False.
     """
-    from rmgpy.molecule.fragment import Fragment
+    if ignore_bond_orders and sanitize:
+        raise ValueError("If ignore_bond_orders is True, sanitize must be False")
+    from rmgpy.molecule.fragment import Fragment, CuttingLabel
     # Sort the atoms before converting to ensure output is consistent
     # between different runs
     if not save_order:
         mol.sort_atoms()
     atoms = mol.vertices
     rd_atom_indices = {}  # dictionary of RDKit atom indices
-    label_dict = {} # store label of atom for Framgent
+    label_dict = {}  # For fragment cutting labels. Key is rdkit atom index, value is label string
     rdkitmol = Chem.rdchem.EditableMol(Chem.rdchem.Mol())
     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
+        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)
         rd_atom.SetNumRadicalElectrons(atom.radical_electrons)
         rd_atom.SetFormalCharge(atom.charge)
-        if atom.element.symbol == 'C' and atom.lone_pairs == 1 and mol.multiplicity == 1: rd_atom.SetNumRadicalElectrons(
-            2)
+        if atom.element.symbol == 'C' and atom.lone_pairs == 1 and mol.multiplicity == 1:
+            rd_atom.SetNumRadicalElectrons(2)
         rdkitmol.AddAtom(rd_atom)
         if remove_h and atom.symbol == 'H':
             pass
         else:
             rd_atom_indices[atom] = index
 
-        # Check if a cutting label is present. If preserve this so that it is added to the SMILES string
-        # Fragment's representative species is Molecule (with CuttingLabel replaced by Si but label as CuttingLabel)
-        # so we use detect_cutting_label to check atom.label
-        _, cutting_label_list = Fragment().detect_cutting_label(atom.label)
-        if cutting_label_list != []:
-            saved_index = index
-            label = atom.label
-            if label in label_dict:
-                label_dict[label].append(saved_index)
-            else:
-                label_dict[label] = [saved_index]
+        # Save cutting labels to add to the SMILES string
+        if atom.label and atom.label in ('R', 'L'):
+            label_dict[index] = atom.label
+
     rd_bonds = Chem.rdchem.BondType
-    orders = {'S': rd_bonds.SINGLE, 'D': rd_bonds.DOUBLE, 'T': rd_bonds.TRIPLE, 'B': rd_bonds.AROMATIC,
-              'Q': rd_bonds.QUADRUPLE}
+    # 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, 'vdW': rd_bonds.ZERO,
+              'H': rd_bonds.HYDROGEN, 'R': rd_bonds.UNSPECIFIED,
+              None: rd_bonds.UNSPECIFIED}
     # Add the bonds
     for atom1 in mol.vertices:
         for atom2, bond in atom1.edges.items():
@@ -106,23 +113,31 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_o
             index1 = atoms.index(atom1)
             index2 = atoms.index(atom2)
             if index1 < index2:
-                order_string = bond.get_order_str()
-                order = orders[order_string]
+                if ignore_bond_orders:
+                    order = rd_bonds.UNSPECIFIED
+                else:
+                    order_string = bond.get_order_str()
+                    order = orders[order_string]
                 rdkitmol.AddBond(index1, index2, order)
 
     # Make editable mol into a mol and rectify the molecule
     rdkitmol = rdkitmol.GetMol()
-    if label_dict:
-        for label, ind_list in label_dict.items():
-            for ind in ind_list:
-                Chem.SetSupplementalSmilesLabel(rdkitmol.GetAtomWithIdx(ind), label)
+    for index, label in label_dict.items():
+        Chem.SetSupplementalSmilesLabel(rdkitmol.GetAtomWithIdx(index), label)
     for atom in rdkitmol.GetAtoms():
         if atom.GetAtomicNum() > 1:
             atom.SetNoImplicit(True)
-    if sanitize:
-        Chem.SanitizeMol(rdkitmol)
     if remove_h:
-        rdkitmol = Chem.RemoveHs(rdkitmol, sanitize=sanitize)
+        rdkitmol = Chem.RemoveHs(rdkitmol, sanitize=False) # skip sanitization here, do it later if requested
+    if sanitize == True:
+        Chem.SanitizeMol(rdkitmol)
+    elif sanitize == "partial":
+        try:
+            Chem.SanitizeMol(rdkitmol, sanitizeOps=Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES)
+        except (KekulizeException, AtomKekulizeException):
+            logging.warning("Kekulization failed; sanitizing without Kekulize")
+            Chem.SanitizeMol(rdkitmol, sanitizeOps=Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES ^ Chem.SANITIZE_KEKULIZE)
+
     if return_mapping:
         return rdkitmol, rd_atom_indices
     return rdkitmol

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,35 @@ 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, saving mapping
+            # in order to match the atoms in the rdmol with the atoms in the
+            # RMG molecule (which is required to extract coordinates).
+            rdmol, rd_atom_idx = self.molecule.to_rdkit_mol(remove_h=False,
+                                                            return_mapping=True,
+                                                            sanitize="partial")
+
+            AllChem.Compute2DCoords(rdmol)
+
+            # Extract the coordinates from each atom.
+            rd_conformer = rdmol.GetConformer(0)
+            for atom in atoms:
+                index = rd_atom_idx[atom]
+                point = rd_conformer.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 +461,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/filtration.py

@@ -401,7 +401,7 @@ def aromaticity_filtration(mol_list, features):
         # Look for structures that don't have standard SDSDSD bond orders
         for mol in other_list:
             # Check all 6 membered rings
-            rings = [ring for ring in mol.get_relevant_cycles() if len(ring) == 6]
+            rings = [ring for ring in mol.get_smallest_set_of_smallest_rings() if len(ring) == 6]
             for ring in rings:
                 bond_list = mol.get_edges_in_cycle(ring)
                 bond_orders = ''.join([bond.get_order_str() for bond in bond_list])