Merge pull request #2329 from ReactionMechanismGenerator/fix_resonanc…

…e_save_order

Fix resonance save order

rmgpy/molecule/converter.pxd

@@ -29,10 +29,10 @@ 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=*)
+cpdef to_rdkit_mol(mm.Molecule mol, bint remove_h=*, bint return_mapping=*, bint sanitize=*, bint save_order=?)
 
 cpdef mm.Molecule from_rdkit_mol(mm.Molecule mol, object rdkitmol, bint raise_atomtype_exception=?)
 
-cpdef to_ob_mol(mm.Molecule mol, bint return_mapping=*)
+cpdef to_ob_mol(mm.Molecule mol, bint return_mapping=*, bint save_order=?)
 
 cpdef mm.Molecule from_ob_mol(mm.Molecule mol, object obmol, bint raise_atomtype_exception=?)

rmgpy/molecule/converter.py

@@ -49,7 +49,7 @@
 from rmgpy.exceptions import DependencyError
 
 
-def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True):
+def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True, save_order=False):
     """
     Convert a molecular structure to a RDKit rdmol object. Uses
     `RDKit <http://rdkit.org/>`_ to perform the conversion.
@@ -61,7 +61,8 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True):
     from rmgpy.molecule.fragment import Fragment
     # Sort the atoms before converting to ensure output is consistent
     # between different runs
-    mol.sort_atoms()
+    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
@@ -227,7 +228,7 @@ def debug_rdkit_mol(rdmol, level=logging.INFO):
     return message
 
 
-def to_ob_mol(mol, return_mapping=False):
+def to_ob_mol(mol, return_mapping=False, save_order=False):
     """
     Convert a molecular structure to an OpenBabel OBMol object. Uses
     `OpenBabel <http://openbabel.org/>`_ to perform the conversion.
@@ -236,7 +237,8 @@ def to_ob_mol(mol, return_mapping=False):
         raise DependencyError('OpenBabel is not installed. Please install or use RDKit.')
 
     # Sort the atoms to ensure consistent output
-    mol.sort_atoms()
+    if not save_order:
+        mol.sort_atoms()
     atoms = mol.vertices
 
     ob_atom_ids = {}  # dictionary of OB atom IDs

rmgpy/molecule/converterTest.py

@@ -117,6 +117,18 @@ def test_atom_mapping_2(self):
                     except RuntimeError:
                         self.fail("RDKit failed in finding the bond in the original atom!")
 
+    def test_atom_mapping_3(self):
+        """Test that to_rdkit_mol with save_order=True retains the atom order and create the correct RDKit Molecule"""
+        adjlist = """1  H u0 p0 c0 {2,S}
+2  C u0 p0 c0 {1,S} {3,T}
+3  N u0 p1 c0 {2,T}
+"""
+        mol = Molecule().from_adjacency_list(adjlist)
+        rdkitmol, _ = to_rdkit_mol(mol, remove_h=False, return_mapping=True, save_order=True)
+
+        self.assertSequenceEqual([atom.number for atom in mol.atoms], [1, 6, 7])
+        self.assertSequenceEqual([rdkitmol.GetAtomWithIdx(idx).GetAtomicNum() for idx in range(3)], [1, 6, 7])
+
 
 class ConverterTest(unittest.TestCase):
 

rmgpy/molecule/fragment.py

@@ -1185,7 +1185,7 @@ def get_representative_molecule(self, mode='minimal', update=True):
 
             return mol_repr, mapping
 
-    def to_rdkit_mol(self, remove_h=False, return_mapping=True):
+    def to_rdkit_mol(self, remove_h=False, return_mapping=True, save_order=False):
         """
         Convert a molecular structure to a RDKit rdmol object.
         """
@@ -1197,17 +1197,19 @@ def to_rdkit_mol(self, remove_h=False, return_mapping=True):
 
         mol0, mapping = self.get_representative_molecule('minimal', update=False)
 
-        rdmol, rdAtomIdx_mol0 = converter.to_rdkit_mol(mol0, remove_h=remove_h,
-                                                     return_mapping=return_mapping,
-                                                     sanitize=True)
+        rdmol, rdAtomIdx_mol0 = converter.to_rdkit_mol(mol0,
+                                                       remove_h=remove_h,
+                                                       return_mapping=return_mapping,
+                                                       sanitize=True,
+                                                       save_order=save_order)
 
         rdAtomIdx_frag = {}
         for frag_atom, mol0_atom in mapping.items():
             rd_idx = rdAtomIdx_mol0[mol0_atom]
             rdAtomIdx_frag[frag_atom] = rd_idx
 
         # sync the order of fragment vertices with the order
-        # of mol0.atoms since mol0.atoms is changed/sorted in 
+        # of mol0.atoms since mol0.atoms is changed/sorted in
         # converter.to_rdkit_mol().
         # Since the rdmol's atoms order is same as the order of mol0's atoms,
         # the synchronization between fragment.atoms order and mol0.atoms order
@@ -1268,7 +1270,7 @@ def from_adjacency_list(self, adjlist, saturate_h=False, raise_atomtype_exceptio
                                  'Currently, AFM does not support ion chemistry.\n {0}'.format(adjlist))
         return self
 
-    def get_aromatic_rings(self, rings=None):
+    def get_aromatic_rings(self, rings=None, save_order=False):
         """
         Returns all aromatic rings as a list of atoms and a list of bonds.
 
@@ -1290,7 +1292,7 @@ def get_aromatic_rings(self, rings=None):
             return [], []
 
         try:
-            rdkitmol, rdAtomIndices = self.to_rdkit_mol(remove_h=False, return_mapping=True)
+            rdkitmol, rdAtomIndices = self.to_rdkit_mol(remove_h=False, return_mapping=True, save_order=save_order)
         except ValueError:
             logging.warning('Unable to check aromaticity by converting to RDKit Mol.')
         else:
@@ -1531,13 +1533,16 @@ def saturate_radicals(self):
 
         return added
 
-    def is_aryl_radical(self, aromatic_rings=None):
+    def is_aryl_radical(self, aromatic_rings=None, save_order=False):
         """
-        Return ``True`` if the fragment only contains aryl radicals,
-        ie. radical on an aromatic ring, or ``False`` otherwise.
+        Return ``True`` if the molecule only contains aryl radicals,
+        i.e., radical on an aromatic ring, or ``False`` otherwise.
+        If no ``aromatic_rings`` provided, aromatic rings will be searched in-place,
+        and this process may involve atom order change by default. Set ``save_order`` to
+        ``True`` to force the atom order unchanged.
         """
         if aromatic_rings is None:
-            aromatic_rings = self.get_aromatic_rings()[0]
+            aromatic_rings = self.get_aromatic_rings(save_order=save_order)[0]
 
         total = self.get_radical_count()
         aromatic_atoms = set([atom for atom in itertools.chain.from_iterable(aromatic_rings)])

rmgpy/molecule/molecule.pxd

@@ -242,16 +242,16 @@ cdef class Molecule(Graph):
     cpdef double calculate_cp0(self) except -1
 
     cpdef double calculate_cpinf(self) except -1
-    
+
     cpdef update_atomtypes(self, bint log_species=?, bint raise_exception=?)
-    
+
     cpdef bint is_radical(self) except -2
 
     cpdef bint has_lone_pairs(self) except -2
 
     cpdef bint has_halogen(self) except -2
 
-    cpdef bint is_aryl_radical(self, list aromatic_rings=?) except -2
+    cpdef bint is_aryl_radical(self, list aromatic_rings=?, bint save_order=?) except -2
 
     cpdef float calculate_symmetry_number(self) except -1
 
@@ -267,7 +267,7 @@ cdef class Molecule(Graph):
 
     cpdef int count_aromatic_rings(self)
 
-    cpdef tuple get_aromatic_rings(self, list rings=?)
+    cpdef tuple get_aromatic_rings(self, list rings=?, bint save_order=?)
 
     cpdef list get_deterministic_sssr(self)
 

rmgpy/molecule/molecule.py

@@ -2140,14 +2140,17 @@ def has_halogen(self):
                 return True
         return False
 
-    def is_aryl_radical(self, aromatic_rings=None):
+    def is_aryl_radical(self, aromatic_rings=None, save_order=False):
         """
         Return ``True`` if the molecule only contains aryl radicals,
-        ie. radical on an aromatic ring, or ``False`` otherwise.
+        i.e., radical on an aromatic ring, or ``False`` otherwise.
+        If no ``aromatic_rings`` provided, aromatic rings will be searched in-place,
+        and this process may involve atom order change by default. Set ``save_order`` to
+        ``True`` to force the atom order unchanged.
         """
         cython.declare(atom=Atom, total=int, aromatic_atoms=set, aryl=int)
         if aromatic_rings is None:
-            aromatic_rings = self.get_aromatic_rings()[0]
+            aromatic_rings = self.get_aromatic_rings(save_order=save_order)[0]
 
         total = self.get_radical_count()
         aromatic_atoms = set([atom for atom in itertools.chain.from_iterable(aromatic_rings)])
@@ -2228,7 +2231,7 @@ def saturate_unfilled_valence(self, update=True):
 
     def saturate_radicals(self, raise_atomtype_exception=True):
         """
-        Saturate the molecule by replacing all radicals with bonds to hydrogen atoms.  Changes self molecule object.  
+        Saturate the molecule by replacing all radicals with bonds to hydrogen atoms.  Changes self molecule object.
         """
         cython.declare(added=dict, atom=Atom, i=int, H=Atom, bond=Bond)
         added = {}
@@ -2341,7 +2344,7 @@ def count_aromatic_rings(self):
 
         return count
 
-    def get_aromatic_rings(self, rings=None):
+    def get_aromatic_rings(self, rings=None, save_order=False):
         """
         Returns all aromatic rings as a list of atoms and a list of bonds.
 
@@ -2351,6 +2354,9 @@ def get_aromatic_rings(self, rings=None):
 
         The method currently restricts aromaticity to six-membered carbon-only rings. This is a limitation imposed
         by RMG, and not by RDKit.
+
+        By default, the atom order will be sorted to get consistent results from different runs. The atom order can
+        be saved when dealing with problems that are sensitive to the atom map.
         """
         cython.declare(rd_atom_indices=dict, ob_atom_ids=dict, aromatic_rings=list, aromatic_bonds=list)
         cython.declare(ring0=list, i=cython.int, atom1=Atom, atom2=Atom)
@@ -2365,7 +2371,7 @@ def get_aromatic_rings(self, rings=None):
             return [], []
 
         try:
-            rdkitmol, rd_atom_indices = converter.to_rdkit_mol(self, remove_h=False, return_mapping=True)
+            rdkitmol, rd_atom_indices = converter.to_rdkit_mol(self, remove_h=False, return_mapping=True, save_order=save_order)
         except ValueError:
             logging.warning('Unable to check aromaticity by converting to RDKit Mol.')
         else:
@@ -2394,7 +2400,7 @@ def get_aromatic_rings(self, rings=None):
 
         logging.info('Trying to use OpenBabel to check aromaticity.')
         try:
-            obmol, ob_atom_ids = converter.to_ob_mol(self, return_mapping=True)
+            obmol, ob_atom_ids = converter.to_ob_mol(self, return_mapping=True, save_order=save_order)
         except DependencyError:
             logging.warning('Unable to check aromaticity by converting for OB Mol.')
             return [], []

rmgpy/molecule/moleculeTest.py

@@ -2501,6 +2501,30 @@ def test_aromaticity_perception_furan(self):
         self.assertEqual(len(aromatic_atoms), 0)
         self.assertEqual(len(aromatic_bonds), 0)
 
+    def test_aromaticity_perception_save_order(self):
+        """Test aromaticity perception via get_aromatic_rings for phenyl radical without changing atom order."""
+        mol = Molecule().from_adjacency_list("""multiplicity 2
+1  C u0 p0 c0 {2,S} {3,S} {7,D}
+2  O u1 p2 c0 {1,S}
+3  C u0 p0 c0 {1,S} {4,D} {8,S}
+4  C u0 p0 c0 {3,D} {5,S} {9,S}
+5  C u0 p0 c0 {4,S} {6,D} {10,S}
+6  C u0 p0 c0 {5,D} {7,S} {11,S}
+7  C u0 p0 c0 {1,D} {6,S} {12,S}
+8  H u0 p0 c0 {3,S}
+9  H u0 p0 c0 {4,S}
+10 H u0 p0 c0 {5,S}
+11 H u0 p0 c0 {6,S}
+12 H u0 p0 c0 {7,S}
+"""
+        )
+        aromatic_atoms, aromatic_bonds = mol.get_aromatic_rings(save_order=True)
+        self.assertEqual(len(aromatic_atoms), 1)
+        self.assertEqual(len(aromatic_bonds), 1)
+        # A quick check for non-changed atom order is to check
+        # if the first atom becomes the oxygen atom after calling `get_aromatic_rings`
+        self.assertFalse(mol.atoms[0].is_oxygen())
+
     def test_aryl_radical_true(self):
         """Test aryl radical perception for phenyl radical."""
         mol = Molecule(smiles='[c]1ccccc1')

rmgpy/molecule/resonance.pxd

@@ -30,7 +30,7 @@ from rmgpy.molecule.molecule cimport Atom, Bond, Molecule
 
 cpdef list populate_resonance_algorithms(dict features=?)
 
-cpdef dict analyze_molecule(Graph mol)
+cpdef dict analyze_molecule(Graph mol, bint save_order=?)
 
 cpdef list generate_resonance_structures(Graph mol, bint clar_structures=?, bint keep_isomorphic=?, bint filter_structures=?, bint save_order=?)
 
@@ -52,14 +52,14 @@ cpdef list generate_isomorphic_resonance_structures(Graph mol, bint saturate_h=?
 
 cpdef list generate_optimal_aromatic_resonance_structures(Graph mol, dict features=?, bint save_order=?)
 
-cpdef list generate_aromatic_resonance_structure(Graph mol, list aromatic_bonds=?, bint copy=?)
+cpdef list generate_aromatic_resonance_structure(Graph mol, list aromatic_bonds=?, bint copy=?, bint save_order=?)
 
 cpdef list generate_aryne_resonance_structures(Graph mol)
 
 cpdef list generate_kekule_structure(Graph mol)
 
-cpdef list generate_clar_structures(Graph mol)
+cpdef list generate_clar_structures(Graph mol, bint save_order=?)
 
-cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?)
+cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?, save_order=?)
 
 cpdef list _clar_transformation(Graph mol, list aromatic_ring)