Merge pull request #1179 from ReactionMechanismGenerator/account_for_…

…enantiomers_in_entropy

Account for enantiomers in entropy

documentation/source/users/rmg/thermo.rst

@@ -266,6 +266,14 @@ Chiral molecules belong to point groups that lack a superposable mirror image
 (i.e. point groups lacking :math:`\sigma_h`, :math:`\sigma_d`, :math:`\sigma_v`, and :math:`S_n`
 symmetry elements).
 
+In RMG, chirality is incorportated into the symmetry attribute by dividing the symmetry by
+two which will increase entropy by :math:`+R * ln(2)`. RMG currently checks for each chiral
+center, defined by 4 different groups attached to a carbon, and halves the symmetry 
+for each chiral center.
+
+The effect of cis-trans isomers is currently not accounted for in RMG.
+
+
 References
 ==========
 
@@ -283,4 +291,4 @@ References
 
 .. [RDKit] Landrum, G. (2012). RDKit, http://rdkit.org.
 
-.. [Benson] Benson, S.W. (1965), https://en.wikipedia.org/wiki/Benson_group_increment_theory
+.. [Benson] Benson, S.W. (1965), https://en.wikipedia.org/wiki/Benson_group_increment_theory

rmgpy/molecule/molecule.pxd

@@ -115,7 +115,7 @@ cdef class Bond(Edge):
 cdef class Molecule(Graph):
 
     cdef public bint implicitHydrogens
-    cdef public int symmetryNumber
+    cdef public float symmetryNumber
     cdef public int multiplicity
     cdef public object rdMol
     cdef public int rdMolConfId
@@ -211,7 +211,7 @@ cdef class Molecule(Graph):
 
     cpdef bint isArylRadical(self, list aromaticRings=?) except -2
 
-    cpdef int calculateSymmetryNumber(self) except -1
+    cpdef float calculateSymmetryNumber(self) except -1
 
     cpdef list generateResonanceIsomers(self, bint keepIsomorphic=?)
 

rmgpy/molecule/molecule.py

@@ -658,7 +658,7 @@ class Molecule(Graph):
     ======================= =========== ========================================
     Attribute               Type        Description
     ======================= =========== ========================================
-    `symmetryNumber`        ``int``     The (estimated) external + internal symmetry number of the molecule
+    `symmetryNumber`        ``float``   The (estimated) external + internal symmetry number of the molecule, modified for chirality
     `multiplicity`          ``int``     The multiplicity of this species, multiplicity = 2*total_spin+1
     `props`                 ``dict``    A list of properties describing the state of the molecule.
     `InChI`                 ``str``     A string representation of the molecule in InChI

rmgpy/molecule/symmetry.pxd

@@ -28,12 +28,12 @@ from .molecule cimport Atom, Bond, Molecule
 
 ################################################################################
 
-cpdef int calculateAtomSymmetryNumber(Molecule molecule, Atom atom) except -1
+cpdef float calculateAtomSymmetryNumber(Molecule molecule, Atom atom) except -1
 
-cpdef int calculateBondSymmetryNumber(Molecule molecule, Atom atom1, Atom atom2) except -1
+cpdef float calculateBondSymmetryNumber(Molecule molecule, Atom atom1, Atom atom2) except -1
 
-cpdef int calculateAxisSymmetryNumber(Molecule molecule) except -1
+cpdef float calculateAxisSymmetryNumber(Molecule molecule) except -1
 
-cpdef int calculateCyclicSymmetryNumber(Molecule molecule) except -1
+cpdef float calculateCyclicSymmetryNumber(Molecule molecule) except -1
 
-cpdef int calculateSymmetryNumber(Molecule molecule) except -1
+cpdef float calculateSymmetryNumber(Molecule molecule) except -1

rmgpy/molecule/symmetry.py

@@ -84,7 +84,7 @@ def calculateAtomSymmetryNumber(molecule, atom):
             elif count == [3, 1]: symmetryNumber *= 3
             elif count == [2, 2]: symmetryNumber *= 2
             elif count == [2, 1, 1]: symmetryNumber *= 1
-            elif count == [1, 1, 1, 1]: symmetryNumber *= 1
+            elif count == [1, 1, 1, 1]: symmetryNumber *= 0.5 # found chirality
         elif single == 3:
             # Three single bonds
             if count == [3]: symmetryNumber *= 3

rmgpy/molecule/symmetryTest.py

@@ -76,6 +76,76 @@ def testAtomSymmetryNumberEthane(self):
                 symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom)
         self.assertEqual(symmetryNumber, 9)
 
+    def testAtomSymmetryNumberEthanewithDeuteriumTritium(self):
+        """
+        Test the Molecule.calculateAtomSymmetryNumber() on CC(D)(T)
+
+        This is meant to test whether chirality is accounted for, which
+        should half the symmetry term.
+
+        The total number is 1.5 because the methyl group contributes *3 and
+        the chiral center contributes *0.5
+        """
+        molecule = Molecule().fromAdjacencyList(
+"""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 C u0 p0 c0 {1,S} {6,S} {7,S} {8,S}
+3 D u0 p0 c0 {1,S}
+4 T u0 p0 c0 {1,S}
+5 H u0 p0 c0 {1,S}
+6 H u0 p0 c0 {2,S}
+7 H u0 p0 c0 {2,S}
+8 H u0 p0 c0 {2,S}
+
+""")
+        symmetryNumber = 1
+        for atom in molecule.atoms:
+            if not molecule.isAtomInCycle(atom):
+                symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom)
+        self.assertAlmostEqual(symmetryNumber, 1.5)
+
+    def testAtomSymmetryNumberWithTwoChiralCenters(self):
+        """
+        Test the Molecule.calculateAtomSymmetryNumber() on [CH2]CC([CH2])C(C)C=C
+
+        This is meant to test whether chirality is accounted for, which
+        should half the symmetry term.
+
+        The molecule has one methyl group (*3), two CH2dot groups (*2 each), and
+        two chiral centers (*0.5), leading to a total atom symmetry number of 3
+        """
+        molecule = Molecule().fromAdjacencyList(
+"""
+multiplicity 3
+1  C u1 p0 c0 {2,S} {3,S} {4,S}
+2  H u0 p0 c0 {1,S}
+3  H u0 p0 c0 {1,S}
+4  C u0 p0 c0 {1,S} {5,S} {13,S} {14,S}
+5  C u0 p0 c0 {4,S} {6,S} {9,S} {15,S}
+6  C u1 p0 c0 {5,S} {7,S} {8,S}
+7  H u0 p0 c0 {6,S}
+8  H u0 p0 c0 {6,S}
+9  C u0 p0 c0 {5,S} {10,S} {11,S} {16,S}
+10 C u0 p0 c0 {9,S} {17,S} {18,S} {19,S}
+11 C u0 p0 c0 {9,S} {12,D} {20,S}
+12 C u0 p0 c0 {11,D} {21,S} {22,S}
+13 H u0 p0 c0 {4,S}
+14 H u0 p0 c0 {4,S}
+15 H u0 p0 c0 {5,S}
+16 H u0 p0 c0 {9,S}
+17 H u0 p0 c0 {10,S}
+18 H u0 p0 c0 {10,S}
+19 H u0 p0 c0 {10,S}
+20 H u0 p0 c0 {11,S}
+21 H u0 p0 c0 {12,S}
+22 H u0 p0 c0 {12,S}
+""")
+        symmetryNumber = 1
+        for atom in molecule.atoms:
+            if not molecule.isAtomInCycle(atom):
+                symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom)
+        self.assertAlmostEqual(symmetryNumber, 3)
+
     def testAtomSymmetryNumberPropane(self):
         """
         Test the Molecule.calculateAtomSymmetryNumber() on CCC

rmgpy/species.pxd

@@ -48,7 +48,7 @@ cdef class Species:
     cdef public object energyTransferModel
     cdef public dict props
     cdef public str aug_inchi
-    cdef public short symmetryNumber
+    cdef public float symmetryNumber
 
     cpdef generateResonanceIsomers(self,bint keepIsomorphic=?)