Merge pull request #1187 from ReactionMechanismGenerator/carbene_constraint

Add carbene species constraint

Author: KEHANG

documentation/source/users/rmg/input.rst

@@ -586,13 +586,15 @@ all of RMG's reaction families. ::
     generatedSpeciesConstraints(
         allowed=['input species','seed mechanisms','reaction libraries'],
         maximumCarbonAtoms=10,
-        maximumHydrogenAtoms=10,
-        maximumOxygenAtoms=10,
-        maximumNitrogenAtoms=10,
-        maximumSiliconAtoms=10,
-        maximumSulfurAtoms=10,
+        maximumOxygenAtoms=2,
+        maximumNitrogenAtoms=2,
+        maximumSiliconAtoms=2,
+        maximumSulfurAtoms=2,
         maximumHeavyAtoms=10,
-        maximumRadicalElectrons=10,
+        maximumRadicalElectrons=2,
+        maximumSingletCarbenes=1,
+        maximumCarbeneRadicals=0,
+        maximumIsotopicAtoms=2,
         allowSingletO2 = False,
     )
 

examples/rmg/commented/input.py

@@ -93,6 +93,16 @@ def failsSpeciesConstraints(species):
         if (struct.getRadicalCount() > maxRadicals):
             return True
 
+    maxCarbenes = speciesConstraints.get('maximumSingletCarbenes', 1)
+    if maxRadicals != -1:
+        if struct.getSingletCarbeneCount() > maxCarbenes:
+            return True
+
+    maxCarbeneRadicals = speciesConstraints.get('maximumCarbeneRadicals', 0)
+    if maxCarbeneRadicals != -1:
+        if struct.getSingletCarbeneCount() > 0 and struct.getRadicalCount() > maxCarbeneRadicals:
+            return True
+
     maxIsotopes = speciesConstraints.get('maximumIsotopicAtoms', -1)
     if maxIsotopes != -1:
         counter = 0

rmgpy/constraints.py

@@ -33,23 +33,210 @@
 """
 
 import unittest
+import mock
 
-from rmgpy.constraints import *
+from rmgpy.rmg.main import RMG
+from rmgpy.constraints import failsSpeciesConstraints
+from rmgpy.species import Species
+from rmgpy.molecule import Molecule
+import rmgpy.rmg.input
 
 ################################################################################
 
 class TestFailsSpeciesConstraints(unittest.TestCase):
     """
     Contains unit tests of the failsSpeciesConstraints function.
     """
-            
-    def testConstraintsNotLoaded(self):
+
+    @classmethod
+    def setUpClass(cls):
+        """
+        A function run ONCE before all unit tests in this class.
+        """
+        cls.rmg = RMG()
+        rmgpy.rmg.input.rmg = cls.rmg
+        rmgpy.rmg.input.generatedSpeciesConstraints(
+            maximumCarbonAtoms=2,
+            maximumOxygenAtoms=1,
+            maximumNitrogenAtoms=1,
+            maximumSiliconAtoms=1,
+            maximumSulfurAtoms=1,
+            maximumHeavyAtoms=3,
+            maximumRadicalElectrons=2,
+            maximumSingletCarbenes=1,
+            maximumCarbeneRadicals=0,
+            maximumIsotopicAtoms=2,
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        """
+        A function run ONCE after all unit tests in this class.
+        """
+        rmgpy.rmg.input.rmg = None
+
+    @mock.patch('rmgpy.constraints.logging')
+    def testConstraintsNotLoaded(self, mock_logging):
         """
         Test what happens when constraints are not loaded.
         """
-        from rmgpy.species import Species
+        # Reset module level rmg variable in rmgpy.rmg.input
+        rmgpy.rmg.input.rmg = None
+
+        mol = Molecule(SMILES='C')
 
+        self.assertFalse(failsSpeciesConstraints(mol))
+
+        mock_logging.debug.assert_called_with('Species constraints could not be found.')
+
+        # Restore module level rmg variable in rmgpy.rmg.input
+        rmgpy.rmg.input.rmg = self.rmg
+
+    def testSpeciesInput(self):
+        """
+        Test that failsSpeciesConstraints can handle a Species object.
+        """
         spc = Species().fromSMILES('C')
 
-        fails = failsSpeciesConstraints(spc)
-        self.assertFalse(fails)
+        self.assertFalse(failsSpeciesConstraints(spc))
+
+    def testExplicitlyAllowedMolecules(self):
+        """
+        Test that we can explicitly allow molecules in species constraints.
+        """
+        mol = Molecule(SMILES='CCCC')
+        self.assertTrue(failsSpeciesConstraints(mol))
+
+        self.rmg.speciesConstraints['explicitlyAllowedMolecules'] = [Molecule(SMILES='CCCC')]
+        self.assertFalse(failsSpeciesConstraints(mol))
+
+    def testCarbonConstraint(self):
+        """
+        Test that we can constrain the max number of carbon atoms.
+        """
+        mol1 = Molecule(SMILES='CC')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='CCC')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testOxygenConstraint(self):
+        """
+        Test that we can constrain the max number of oxygen atoms.
+        """
+        mol1 = Molecule(SMILES='C=O')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='OC=O')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testNitrogenConstraint(self):
+        """
+        Test that we can constrain the max number of nitrogen atoms.
+        """
+        mol1 = Molecule(SMILES='CN')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='NCN')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testSiliconConstraint(self):
+        """
+        Test that we can constrain the max number of silicon atoms.
+        """
+        mol1 = Molecule(SMILES='[SiH4]')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='[SiH3][SiH3]')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testSulfurConstraint(self):
+        """
+        Test that we can constrain the max number of sulfur atoms.
+        """
+        mol1 = Molecule(SMILES='CS')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='SCS')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testHeavyConstraint(self):
+        """
+        Test that we can constrain the max number of heavy atoms.
+        """
+        mol1 = Molecule(SMILES='CCO')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='CCN=O')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testRadicalConstraint(self):
+        """
+        Test that we can constrain the max number of radical electrons.
+        """
+        mol1 = Molecule(SMILES='[CH2][CH2]')
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule(SMILES='[CH2][CH][CH2]')
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testCarbeneConstraint(self):
+        """
+        Test that we can constrain the max number of singlet carbenes.
+        """
+        mol1 = Molecule().fromAdjacencyList("""
+1 C u0 p1 c0 {2,S} {3,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+""")
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule().fromAdjacencyList("""
+1 C u0 p1 c0 {2,S} {3,S}
+2 H u0 p0 c0 {1,S}
+3 C u0 p1 c0 {1,S} {4,S}
+4 H u0 p0 c0 {3,S}
+""")
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testCarbeneRadicalConstraint(self):
+        """
+        Test that we can constrain the max number of radical electrons with a carbene.
+        """
+        mol1 = Molecule().fromAdjacencyList("""
+1 C u0 p1 c0 {2,S} {3,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+""")
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule().fromAdjacencyList("""
+1 C u0 p1 c0 {2,S} {3,S}
+2 H u0 p0 c0 {1,S}
+3 C u1 p0 c0 {1,S} {4,S} {5,S}
+4 H u0 p0 c0 {3,S}
+5 H u0 p0 c0 {3,S}
+""")
+        self.assertTrue(failsSpeciesConstraints(mol2))
+
+    def testIsotopeConstraint(self):
+        """
+        Test that we can constrain the max number of isotopic atoms.
+        """
+        mol1 = Molecule().fromAdjacencyList("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 D u0 p0 c0 {1,S}
+3 D u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 H u0 p0 c0 {1,S}
+""")
+        self.assertFalse(failsSpeciesConstraints(mol1))
+
+        mol2 = Molecule().fromAdjacencyList("""
+1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 D u0 p0 c0 {1,S}
+3 D u0 p0 c0 {1,S}
+4 D u0 p0 c0 {1,S}
+5 H u0 p0 c0 {1,S}
+""")
+        self.assertTrue(failsSpeciesConstraints(mol2))

rmgpy/constraintsTest.py

@@ -1329,8 +1329,11 @@ def isMoleculeForbidden(self, molecule):
 
         forbidden_structures = getDB('forbidden')
 
+        # check family-specific forbidden structures 
         if self.forbidden is not None and self.forbidden.isMoleculeForbidden(molecule):
             return True
+
+        # check RMG globally forbidden structures
         if forbidden_structures.isMoleculeForbidden(molecule):
             return True
         return False

rmgpy/data/kinetics/family.py

@@ -145,6 +145,8 @@ cdef class Molecule(Graph):
 
     cpdef short getRadicalCount(self)
 
+    cpdef short getSingletCarbeneCount(self)
+
     cpdef double getMolecularWeight(self)
 
     cpdef int getNumAtoms(self, str element=?)

rmgpy/molecule/molecule.pxd

@@ -884,6 +884,19 @@ def getRadicalCount(self):
             radicals += atom.radicalElectrons
         return radicals
 
+    def getSingletCarbeneCount(self):
+        """
+        Return the total number of singlet carbenes (lone pair on a carbon atom)
+        in the molecule. Counts the number of carbon atoms with a lone pair.
+        In the case of [C] with two lone pairs, this method will return 1.
+        """
+        cython.declare(atom=Atom, carbenes=cython.short)
+        carbenes = 0
+        for atom in self.vertices:
+            if atom.isCarbon() and atom.lonePairs > 0:
+                carbenes += 1
+        return carbenes
+
     def getNumAtoms(self, element = None):
         """
         Return the number of atoms in molecule.  If element is given, ie. "H" or "C",

rmgpy/molecule/molecule.py

@@ -1089,6 +1089,29 @@ def testRadicalCH2CH2CH2(self):
         """
         molecule = Molecule().fromSMILES('[CH2]C[CH2]')
         self.assertEqual(molecule.getRadicalCount(), 2)
+
+    def testSingletCarbene(self):
+        """Test radical and carbene count on singlet carbene."""
+        mol = Molecule().fromAdjacencyList("""
+1 C u0 p1 {2,S}
+2 C u0 p1 {1,S}
+""", saturateH=True)
+        self.assertEqual(mol.getRadicalCount(), 0)
+        self.assertEqual(mol.getSingletCarbeneCount(), 2)
+
+    def testTripletCarbene(self):
+        """Test radical and carbene count on triplet carbene."""
+        mol = Molecule().fromAdjacencyList("""
+1 C u2 p0 {2,S}
+2 C u0 p1 {1,S}
+""", saturateH=True)
+        self.assertEqual(mol.getRadicalCount(), 2)
+        self.assertEqual(mol.getSingletCarbeneCount(), 1)
+
+    def testSingletCarbon(self):
+        """Test that getSingletCarbeneCount returns 1 for singlet carbon atom."""
+        mol = Molecule().fromAdjacencyList('1 C u0 p2')
+        self.assertEqual(mol.getSingletCarbeneCount(), 1)
 
     def testSMILES(self):
         """

rmgpy/molecule/moleculeTest.py

@@ -353,6 +353,8 @@ def generatedSpeciesConstraints(**kwargs):
         'maximumSulfurAtoms',
         'maximumHeavyAtoms',
         'maximumRadicalElectrons',
+        'maximumSingletCarbenes',
+        'maximumCarbeneRadicals',
         'allowSingletO2',
         'maximumIsotopicAtoms'
     ]