From 3a73cff8d153e273158ea9f53e3e3cff90c9e71f Mon Sep 17 00:00:00 2001
From: Enoch Dames <enoch.dames@gmail.com>
Date: Fri, 16 May 2014 19:05:00 -0400
Subject: [PATCH] added aromaticity recognition for symmetry calcs

Aromaticity is now perceived in calculating cyclic symmetry numbers.
This commit copies a molecule inside `calculateCyclicSymmetryNumber`,
turns it into an rdkit object, finds all the aromatic bonds and atoms.
The same bonds and atoms in the corresponding rmg molecule object are
then redifined so that the alogorithm is using the correct information.

There is one strange thing about this fix that I dont yet understand.
The unittest thinks cyclic symmetry number of dimethylbenzene is 1
(should be 2). On my local machine, the same
`calculateCyclicSymmetryNumber` function gets it correct, at 2. Also,
i'm inadvertantly printing a bit to stdout, and not sure why it's
happening?

relevant to issues #12, #24, #119, #70, #141
---
 rmgpy/molecule/molecule.py     |  2 ++
 rmgpy/molecule/symmetry.py     | 20 +++++++++++++++++---
 rmgpy/molecule/symmetryTest.py |  1 -
 3 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index cd1cac7104..9e69e06a38 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -567,6 +567,8 @@ def applyAction(self, action):
                 self.incrementOrder()
             elif action[2] == -1:
                 self.decrementOrder()
+            elif action[2] == 'B':
+                self.order = 'B'
             else:
                 raise ActionError('Unable to update Bond due to CHANGE_BOND action: Invalid order "{0}".'.format(action[2]))
         else:
diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py
index c9dd35e0f2..e938e06dc9 100644
--- a/rmgpy/molecule/symmetry.py
+++ b/rmgpy/molecule/symmetry.py
@@ -343,9 +343,12 @@ def calculateCyclicSymmetryNumber(molecule):
     Get the symmetry number correction for cyclic regions of a molecule.
     For complicated fused rings the smallest set of smallest rings is used.
     """
-
+    from rdkit.Chem.rdmolops import SanitizeMol
+    from rdkit.Chem.rdchem import Mol 
+    mcopy = molecule.toRDKitMol(removeHs=True, returnMapping=False)
+    SanitizeMol(mcopy)
     symmetryNumber = 1
-
+            
     # Get symmetry number for each ring in structure
     rings = molecule.getSmallestSetOfSmallestRings()
     for ring0 in rings:
@@ -353,7 +356,18 @@ def calculateCyclicSymmetryNumber(molecule):
         # Make copy of structure
         structure = molecule.copy(True)
         ring = [structure.atoms[molecule.atoms.index(atom)] for atom in ring0]
-
+        # Figure out which atoms and bonds are aromatic and reassign appropriately:
+        for i, atom1 in enumerate(ring0):
+            for atom2 in ring0[i+1:]:
+                if molecule.hasBond(atom1, atom2):
+                    if str(mcopy.GetBondBetweenAtoms(i,i+1)) != 'None':
+                        if str(mcopy.GetBondBetweenAtoms(i,i+1).GetBondType()) == 'AROMATIC':
+                            bond = molecule.getBond(atom1, atom2)
+                            bond.applyAction(['CHANGE_BOND', atom1, 'B', atom2])
+                            atom1.atomType.label = 'Cb'
+                            atom2.atomType.label = 'Cb'
+                    else:
+                        pass
         # Remove bonds of ring from structure
         for i, atom1 in enumerate(ring):
             for atom2 in ring[i+1:]:
diff --git a/rmgpy/molecule/symmetryTest.py b/rmgpy/molecule/symmetryTest.py
index 953a857d9c..533f3750d4 100644
--- a/rmgpy/molecule/symmetryTest.py
+++ b/rmgpy/molecule/symmetryTest.py
@@ -244,7 +244,6 @@ def testCyclicSymmetryNumberCyclohexane(self):
         symmetryNumber = calculateCyclicSymmetryNumber(molecule)
         self.assertEqual(symmetryNumber, 12)
 
-    @work_in_progress
     def testCyclicSymmetryNumberBenzene(self):
         """
         Test the Molecule.calculateCyclicSymmetryNumber() method.