Merge pull request #208 from connie/master

Modifying species constraints checks

Author: rwest

documentation/source/users/rmg/database/kinetics.rst

@@ -162,6 +162,46 @@ Additionally, groups can also be defined as unions of other groups. For example,
 
 	label="X_H_or_Xrad_H",
 	group=OR{X_H, Xrad_H}, 
+    
+
+Forbidden Groups
+----------------
+Forbidden groups can be defined to ban structures globally in RMG or to
+ban pathways in a specific kinetic family.
+
+Globally forbidden structures will ban all reactions containing either reactants
+or products that are forbidden.  These groups are stored in in the file located at
+``RMG-database/input/forbiddenStructures.py``. 
+
+
+To ban certain specific pathways in the kinetics 
+families, a `forbidden` group must be created, like the following group
+in the ``intra_H_migration`` family ::
+
+    forbidden(
+        label = "bridged56_1254",
+    group =
+    """""""
+    1 *1 C 1 {2,S} {6,S}
+    2 *4 C 0 {1,S} {3,S} {7,S}
+    3    C 0 {2,S} {4,S}
+    4 *2 C 0 {3,S} {5,S} {8,S}
+    5 *5 C 0 {4,S} {6,S} {7,S}
+    6    C 0 {1,S} {5,S}
+    7    C 0 {2,S} {5,S}
+    8 *3 H 0 {4,S}
+    """,
+        shortDesc = u"""""",
+        longDesc = 
+    u"""
+    
+    """,
+    )
+
+Forbidden groups should be placed inside the groups.py file located inside the
+specific kinetics family's folder ``RMG-database/input/kinetics/family_name/`` 
+alongside normal group entries. The starred atoms in the forbidden group
+ban the specified reaction recipe from occurring in matched products and reactants.
 
 Hierarchical Trees
 ------------------

documentation/source/users/rmg/input.rst

@@ -360,14 +360,22 @@ Miscellaneous Options
 
 Miscellaneous options:: 
 
-	options(
-		units='si',
-		saveRestartPeriod=(1,'hour'),
-		drawMolecules=False,
-		generatePlots=False,
-	)
+    options(
+        units='si',
+        saveRestartPeriod=(1,'hour'),
+        drawMolecules=False,
+        generatePlots=False,
+    )
+    
+Species Constraints
+===================== 
 
-	generatedSpeciesConstraints(
+RMG can generate mechanisms with a number of optional species constraints,
+such as total number of carbon atoms or electrons per species. These are applied to
+all of RMG's reaction families. ::
+
+    generatedSpeciesConstraints(
+        allowed=['input species','seed mechanisms','reaction libraries'],
         maximumCarbonAtoms=10,
         maximumHydrogenAtoms=10,
         maximumOxygenAtoms=10,
@@ -376,7 +384,18 @@ Miscellaneous options::
         maximumSulfurAtoms=10,
         maximumHeavyAtoms=10,
         maximumRadicalElectrons=10,
-	)
+    )
+
+An additional flag ``allowed`` can be set to allow species 
+from either the input file, seed mechanisms, or reaction libraries to bypass these constraints.
+Note that this should be done with caution, since the constraints will still apply to subsequent
+products that form.  
+
+Note that under all circumstances all forbidden species will still be banned unless they are 
+manually removed from the database.  See :ref:`kineticsDatabase` for more information on 
+forbidden groups.  
+
+
 
 Examples
 ========

rmgpy/data/kinetics/__init__.py

@@ -62,7 +62,8 @@ def __init__(self):
         self.recommendedFamilies = {}
         self.families = {}
         self.libraries = {}
-        self.libraryOrder = []
+        self.libraryOrder = []     # a list of tuples in the format ('library_label', LibraryType),
+                                   # where LibraryType is set to either 'Reaction Library' or 'Seed'.  
         self.local_context = {
             'KineticsData': KineticsData,
             'Arrhenius': Arrhenius,
@@ -201,7 +202,7 @@ def loadLibraries(self, path, libraries=None):
         The `path` points to the folder of kinetics libraries in the database,
         and the libraries should be in files like :file:`<path>/<library>.py`.
         """
-        self.libraries = {}; self.libraryOrder = []
+        self.libraries = {}
 
         if libraries is not None:
             for library_name in libraries:
@@ -211,11 +212,12 @@ def loadLibraries(self, path, libraries=None):
                     library = KineticsLibrary(label=library_name)
                     library.load(library_file, self.local_context, self.global_context)
                     self.libraries[library.label] = library
-                    self.libraryOrder.append(library.label)
                 else:
                     raise IOError("Couldn't find kinetics library {0}".format(library_file))
-            assert (self.libraryOrder == libraries)
-        else:# load all the libraries you can find
+            # library order should've been set prior to this, with the given seed mechs and reaction libraries
+            assert (len(self.libraryOrder) == len(libraries))
+        else:# load all the libraries you can find (this cannot be activated in a normal RMG job.  Only activated when loading the database for other purposes)
+            self.libraryOrder = []
             for (root, dirs, files) in os.walk(os.path.join(path)):
                 for f in files:
                     name, ext = os.path.splitext(f)
@@ -226,7 +228,7 @@ def loadLibraries(self, path, libraries=None):
                         library = KineticsLibrary(label=label)
                         library.load(library_file, self.local_context, self.global_context)
                         self.libraries[library.label] = library
-                        self.libraryOrder.append(library.label)
+                        self.libraryOrder.append((library.label,'Reaction Library'))
 
     def save(self, path):
         """
@@ -330,29 +332,31 @@ def saveOld(self, path):
                 onoff = 'on ' if self.recommendedFamilies[label] else 'off'
                 f.write("{num:<2d}    {onoff}     {label}\n".format(num=number, label=label, onoff=onoff))
 
-    def generateReactions(self, reactants, products=None, **options):
+    def generateReactions(self, reactants, products=None, failsSpeciesConstraints=None):
         """
         Generate all reactions between the provided list of one or two
         `reactants`, which should be :class:`Molecule` objects. This method
         searches the depository, libraries, and groups, in that order.
         """
         reactionList = []
-        reactionList.extend(self.generateReactionsFromLibraries(reactants, products, **options))
-        reactionList.extend(self.generateReactionsFromFamilies(reactants, products, **options))
+        reactionList.extend(self.generateReactionsFromLibraries(reactants, products, failsSpeciesConstraints=failsSpeciesConstraints))
+        reactionList.extend(self.generateReactionsFromFamilies(reactants, products, failsSpeciesConstraints=failsSpeciesConstraints))
         return reactionList
 
-    def generateReactionsFromLibraries(self, reactants, products, **options):
+    def generateReactionsFromLibraries(self, reactants, products, failsSpeciesConstraints=None):
         """
         Generate all reactions between the provided list of one or two
         `reactants`, which should be :class:`Molecule` objects. This method
         searches the depository.
         """
         reactionList = []
-        for label in self.libraryOrder:
-            reactionList.extend(self.generateReactionsFromLibrary(reactants, products, self.libraries[label], **options))
+        for label, libraryType in self.libraryOrder:
+            # Generate reactions from reaction libraries (no need to generate them from seeds)
+            if libraryType == "Reaction Library":
+                reactionList.extend(self.generateReactionsFromLibrary(reactants, products, self.libraries[label], failsSpeciesConstraints=failsSpeciesConstraints))
         return reactionList
 
-    def generateReactionsFromLibrary(self, reactants, products, library, **options):
+    def generateReactionsFromLibrary(self, reactants, products, library, failsSpeciesConstraints=None):
         """
         Generate all reactions between the provided list of one or two
         `reactants`, which should be :class:`Molecule` objects. This method
@@ -375,7 +379,7 @@ def generateReactionsFromLibrary(self, reactants, products, library, **options):
             reactionList = filterReactions(reactants, products, reactionList)
         return reactionList
 
-    def generateReactionsFromFamilies(self, reactants, products, only_families=None, **options):
+    def generateReactionsFromFamilies(self, reactants, products, only_families=None, failsSpeciesConstraints=None):
         """
         Generate all reactions between the provided list of one or two
         `reactants`, which should be :class:`Molecule` objects. This method
@@ -392,7 +396,7 @@ def generateReactionsFromFamilies(self, reactants, products, only_families=None,
         reactionList = []
         for label, family in self.families.iteritems():
             if only_families is None or label in only_families:
-                reactionList.extend(family.generateReactions(reactants, **options))
+                reactionList.extend(family.generateReactions(reactants, failsSpeciesConstraints=failsSpeciesConstraints))
         if products:
             reactionList = filterReactions(reactants, products, reactionList)
         return reactionList

rmgpy/data/kinetics/family.py

@@ -1078,7 +1078,7 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True):
         # Return the product structures
         return productStructures
 
-    def __generateProductStructures(self, reactantStructures, maps, forward, **options):
+    def __generateProductStructures(self, reactantStructures, maps, forward, failsSpeciesConstraints=None):
         """
         For a given set of `reactantStructures` and a given set of `maps`,
         generate and return the corresponding product structures. The
@@ -1087,6 +1087,8 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio
         parameter is a list of mappings of the top-level tree node of each
         *template* reactant to the corresponding *structure*. This function
         returns a list of the product structures.
+        `failsSpeciesConstraints` is a function that accepts a :class:`Molecule`
+        structure and returns `True` if it is forbidden.
         """
 
         productStructuresList = []
@@ -1125,33 +1127,11 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio
                 productStructures.reverse()
 
         # Apply the generated species constraints (if given)
-        if options:
-            maxCarbonAtoms = options.get('maximumCarbonAtoms', 1000000)
-            maxHydrogenAtoms = options.get('maximumHydrogenAtoms', 1000000)
-            maxOxygenAtoms = options.get('maximumOxygenAtoms', 1000000)
-            maxNitrogenAtoms = options.get('maximumNitrogenAtoms', 1000000)
-            maxSiliconAtoms = options.get('maximumSiliconAtoms', 1000000)
-            maxSulfurAtoms = options.get('maximumSulfurAtoms', 1000000)
-            maxHeavyAtoms = options.get('maximumHeavyAtoms', 1000000)
-            maxRadicals = options.get('maximumRadicalElectrons', 1000000)
+        if failsSpeciesConstraints:
             for struct in productStructures:
-                H = struct.getNumAtoms('H')
-                if struct.getNumAtoms('C') > maxCarbonAtoms:
-                    raise ForbiddenStructureException()
-                if H > maxHydrogenAtoms:
-                    raise ForbiddenStructureException()
-                if struct.getNumAtoms('O') > maxOxygenAtoms:
-                    raise ForbiddenStructureException()
-                if struct.getNumAtoms('N') > maxNitrogenAtoms:
-                    raise ForbiddenStructureException()
-                if struct.getNumAtoms('Si') > maxSiliconAtoms:
-                    raise ForbiddenStructureException()
-                if struct.getNumAtoms('S') > maxSulfurAtoms:
-                    raise ForbiddenStructureException()
-                if len(struct.atoms) - H > maxHeavyAtoms:
-                    raise ForbiddenStructureException()
-                if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H > 1):
-                    raise ForbiddenStructureException()
+                if failsSpeciesConstraints(struct):
+                    raise ForbiddenStructureException() 
+
 
         # Generate other possible electronic states
         electronicStructuresList1 = []
@@ -1389,7 +1369,7 @@ def __matchReactantToTemplate(self, reactant, templateReactant):
         elif isinstance(struct, Group):
             return reactant.findSubgraphIsomorphisms(struct)
 
-    def generateReactions(self, reactants, **options):
+    def generateReactions(self, reactants, failsSpeciesConstraints=None):
         """
         Generate all reactions between the provided list of one or two
         `reactants`, which should be either single :class:`Molecule` objects
@@ -1402,12 +1382,12 @@ def generateReactions(self, reactants, **options):
         reactionList = []
 
         # Forward direction (the direction in which kinetics is defined)
-        reactionList.extend(self.__generateReactions(reactants, forward=True, **options))
+        reactionList.extend(self.__generateReactions(reactants, forward=True, failsSpeciesConstraints=failsSpeciesConstraints))
 
         if self.ownReverse:
             # for each reaction, make its reverse reaction and store in a 'reverse' attribute
             for rxn in reactionList:
-                reactions = self.__generateReactions(rxn.products, products=rxn.reactants, forward=True, **options)
+                reactions = self.__generateReactions(rxn.products, products=rxn.reactants, forward=True, failsSpeciesConstraints=failsSpeciesConstraints)
                 if len(reactions) != 1:
                     logging.error("Expecting one matching reverse reaction, not {0} in reaction family {1} for forward reaction {2}.\n".format(len(reactions), self.label, str(rxn)))
                     for reactant in rxn.reactants:
@@ -1421,7 +1401,7 @@ def generateReactions(self, reactants, **options):
 
         else: # family is not ownReverse
             # Reverse direction (the direction in which kinetics is not defined)
-            reactionList.extend(self.__generateReactions(reactants, forward=False, **options))
+            reactionList.extend(self.__generateReactions(reactants, forward=False, failsSpeciesConstraints=failsSpeciesConstraints))
 
         # Return the reactions as containing Species objects, not Molecule objects
         for reaction in reactionList:
@@ -1450,14 +1430,16 @@ def calculateDegeneracy(self, reaction):
             raise Exception('Unable to calculate degeneracy for reaction {0} in reaction family {1}.'.format(reaction, self.label))
         return reactions[0].degeneracy
 
-    def __generateReactions(self, reactants, products=None, forward=True, **options):
+    def __generateReactions(self, reactants, products=None, forward=True, failsSpeciesConstraints=None):
         """
         Generate a list of all of the possible reactions of this family between
         the list of `reactants`. The number of reactants provided must match
         the number of reactants expected by the template, or this function
         will return an empty list. Each item in the list of reactants should
         be a list of :class:`Molecule` objects, each representing a resonance
         isomer of the species of interest.
+        `failsSpeciesConstraints` is an optional function that accepts a :class:`Molecule`
+        structure and returns `True` if it is forbidden.
         """
 
         rxnList = []; speciesList = []
@@ -1493,7 +1475,7 @@ def __generateReactions(self, reactants, products=None, forward=True, **options)
                 for map in mappings:
                     reactantStructures = [molecule]
                     try:
-                        productStructuresList = self.__generateProductStructures(reactantStructures, [map], forward, **options)
+                        productStructuresList = self.__generateProductStructures(reactantStructures, [map], forward, failsSpeciesConstraints=failsSpeciesConstraints)
                     except ForbiddenStructureException:
                         pass
                     else:
@@ -1521,7 +1503,7 @@ def __generateReactions(self, reactants, products=None, forward=True, **options)
                         for mapB in mappingsB:
                             reactantStructures = [moleculeA, moleculeB]
                             try:
-                                productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options)
+                                productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, failsSpeciesConstraints=failsSpeciesConstraints)
                             except ForbiddenStructureException:
                                 pass
                             else:
@@ -1542,7 +1524,7 @@ def __generateReactions(self, reactants, products=None, forward=True, **options)
                             for mapB in mappingsB:
                                 reactantStructures = [moleculeA, moleculeB]
                                 try:
-                                    productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options)
+                                    productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, failsSpeciesConstraints=failsSpeciesConstraints)
                                 except ForbiddenStructureException:
                                     pass
                                 else:

rmgpy/data/rmg.py

@@ -139,17 +139,20 @@ def loadKinetics(self,
         `path` points to the top-level folder of the RMG kinetics database.
         """
         kineticsLibraries = []
-        if reactionLibraries is not None and seedMechanisms is not None:
-            kineticsLibraries.extend(seedMechanisms)
-            kineticsLibraries.extend(reactionLibraries)
-        elif reactionLibraries is not None and seedMechanisms is None:
-            kineticsLibraries.extend(reactionLibraries)
-        elif reactionLibraries is None and seedMechanisms is not None:
-            kineticsLibraries.extend(seedMechanisms)
-        else:
+        libraryOrder = []
+        if seedMechanisms is None and reactionLibraries is None:
             kineticsLibraries = None
-            
+        if seedMechanisms is not None:
+            for library in seedMechanisms:
+                kineticsLibraries.append(library)
+                libraryOrder.append((library,'Seed'))
+        if reactionLibraries is not None:
+            for library in reactionLibraries:
+                kineticsLibraries.append(library)
+                libraryOrder.append((library,'Reaction Library'))
+        
         self.kinetics = KineticsDatabase()
+        self.kinetics.libraryOrder = libraryOrder
         self.kinetics.load(path,
                            families=kineticsFamilies,
                            libraries=kineticsLibraries,

rmgpy/rmg/input.py

@@ -253,6 +253,7 @@ def options(units='si', saveRestartPeriod=None, drawMolecules=False, generatePlo
 
 def generatedSpeciesConstraints(**kwargs):
     validConstraints = [
+        'allowed',
         'maximumCarbonAtoms',
         'maximumHydrogenAtoms',
         'maximumOxygenAtoms',
@@ -266,7 +267,7 @@ def generatedSpeciesConstraints(**kwargs):
     for key, value in kwargs.items():
         if key not in validConstraints:
             raise InputError('Invalid generated species constraint {0!r}.'.format(key))
-        rmg.reactionGenerationOptions[key] = value
+        rmg.speciesConstraints[key] = value
 
 ################################################################################