Reducing code duplication by passing in failsSpeciesConstraints.

Now passes species constraints function through to family reaction
generation routines instead of running checks in multiple locations.

Author: connie

rmgpy/data/kinetics/__init__.py

@@ -332,18 +332,18 @@ 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
@@ -353,11 +353,10 @@ def generateReactionsFromLibraries(self, reactants, products, **options):
         for label, libraryType in self.libraryOrder:
             # Generate reactions from reaction libraries (no need to generate them from seeds)
             if libraryType == "Reaction Library":
-                print 'generating reactions from library {0}'.format(label)
-                reactionList.extend(self.generateReactionsFromLibrary(reactants, products, self.libraries[label], **options))
+                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
@@ -380,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
@@ -397,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,41 +1127,11 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio
                 productStructures.reverse()
 
         # Apply the generated species constraints (if given)
-        if options:
-            explicitlyAllowedMolecules = options.get('explicitlyAllowedMolecules')
-            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:
-                allowed = False
-                for molecule in explicitlyAllowedMolecules:
-                    if struct.isIsomorphic(molecule):                        
-                        allowed = True
-                        break
-                if allowed:
-                    continue
-                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 = []
@@ -1397,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
@@ -1410,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:
@@ -1429,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:
@@ -1458,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 = []
@@ -1501,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:
@@ -1529,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:
@@ -1550,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/rmg/model.py

@@ -603,15 +603,14 @@ def react(self, database, speciesA, speciesB=None):
         Generates reactions involving :class:`rmgpy.species.Species` speciesA and speciesB.
         """
         reactionList = []
-        options = self.speciesConstraints
         if speciesB is None:
             for moleculeA in speciesA.molecule:
-                reactionList.extend(database.kinetics.generateReactions([moleculeA], **options))
+                reactionList.extend(database.kinetics.generateReactions([moleculeA], failsSpeciesConstraints=self.failsSpeciesConstraints))
                 moleculeA.clearLabeledAtoms()
         else:
             for moleculeA in speciesA.molecule:
                 for moleculeB in speciesB.molecule:
-                    reactionList.extend(database.kinetics.generateReactions([moleculeA, moleculeB], **options))
+                    reactionList.extend(database.kinetics.generateReactions([moleculeA, moleculeB], failsSpeciesConstraints=self.failsSpeciesConstraints))
                     moleculeA.clearLabeledAtoms()
                     moleculeB.clearLabeledAtoms()
         return reactionList
@@ -1706,10 +1705,10 @@ def saveChemkinFile(self, path, verbose_path, dictionaryPath=None, transportPath
 
     def failsSpeciesConstraints(self, species):
         """
-        Checks whether the species passes the speciesConstraints set by the user.  If not,
-        returns `True` for failing speciesConstraints.
+        Pass in either a `Species` or `Molecule` object and checks whether it passes 
+        the speciesConstraints set by the user.  If not, returns `True` for failing speciesConstraints.
         """
-        explicitlyAllowedMolecules = self.speciesConstraints.get('explicitlyAllowedMolecules')
+        explicitlyAllowedMolecules = self.speciesConstraints.get('explicitlyAllowedMolecules', [])
         maxCarbonAtoms = self.speciesConstraints.get('maximumCarbonAtoms', 1000000)
         maxHydrogenAtoms = self.speciesConstraints.get('maximumHydrogenAtoms', 1000000)
         maxOxygenAtoms = self.speciesConstraints.get('maximumOxygenAtoms', 1000000)
@@ -1719,7 +1718,11 @@ def failsSpeciesConstraints(self, species):
         maxHeavyAtoms = self.speciesConstraints.get('maximumHeavyAtoms', 1000000)
         maxRadicals = self.speciesConstraints.get('maximumRadicalElectrons', 1000000)
 
-        struct = species.molecule[0]
+        if isinstance(species, rmgpy.species.Species):
+            struct = species.molecule[0]
+        else:
+            # expects a molecule here
+            struct = species
         for molecule in explicitlyAllowedMolecules:
             if struct.isIsomorphic(molecule):
                 return False