Add species constraints checks to input species, seed, and reaction libraries.

Explicitly allow species by using the `allowed` flag in generatedSpeciesConstraints, which
can be set to a list with 'input species','reaction libraries', or 'seed mechanisms'
(or any combination of these categories).  The species from the specified categories can bypass and
ignore species constraints given in the input file. They are also added to a list
of explicitlyAllowedMolecules.  However, all family forbidden and globally forbidden
groups are still forbidden.

Additionally, now prints warnings for globally forbidden but allowed molecules.  And
gives exceptions if they were not allowed to begin with.

Author: connie

rmgpy/data/kinetics/family.py

@@ -1126,6 +1126,7 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio
 
         # 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)
@@ -1135,6 +1136,13 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio
             maxHeavyAtoms = options.get('maximumHeavyAtoms', 1000000)
             maxRadicals = options.get('maximumRadicalElectrons', 1000000)
             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()

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',

rmgpy/rmg/main.py

@@ -50,6 +50,7 @@
 from rmgpy.solver.base import TerminationTime, TerminationConversion
 from rmgpy.solver.simple import SimpleReactor
 from rmgpy.data.rmg import RMGDatabase
+from rmgpy.data.base import ForbiddenStructureException
 from rmgpy.data.kinetics import KineticsLibrary, KineticsFamily, LibraryReaction, TemplateReaction
 
 from rmgpy.reaction import Reaction
@@ -174,6 +175,7 @@ def loadInput(self, path=None):
         from input import readInputFile
         if path is None: path = self.inputFile
         readInputFile(path, self)
+        self.speciesConstraints['explicitlyAllowedMolecules'] = [] 
         self.reactionModel.kineticsEstimator = self.kineticsEstimator
         # If the output directory is not yet set, then set it to the same
         # directory as the input file by default
@@ -346,6 +348,20 @@ def initialize(self, args):
             for library, option in self.reactionLibraries:
                 self.reactionModel.addReactionLibraryToEdge(library)
 
+            # Perform species constraints and forbidden species checks on input species
+            for spec in self.initialSpecies:
+                if self.database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
+                    if 'allowed' in self.speciesConstraints and 'input species' in self.speciesConstraints['allowed']:
+                        logging.warning('Input species {0} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.'.format(spec.label))
+                    else:
+                        raise ForbiddenStructureException("Input species {0} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label))
+                if self.reactionModel.failsSpeciesConstraints(spec):
+                    if 'allowed' in self.speciesConstraints and 'input species' in self.speciesConstraints['allowed']:
+                        self.speciesConstraints['explicitlyAllowedMolecules'].append(spec.molecule[0])
+                        pass
+                    else:
+                        raise ForbiddenStructureException("Species constraints forbids input species {0}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label))
+            
             # Add nonreactive species (e.g. bath gases) to core first
             # This is necessary so that the PDep algorithm can identify the bath gas
             self.reactionModel.enlarge([spec for spec in self.initialSpecies if not spec.reactive])

rmgpy/rmg/model.py

@@ -47,7 +47,7 @@
 from rmgpy.pdep import SingleExponentialDown
 from rmgpy.statmech import  Conformer
 
-from rmgpy.data.base import Entry
+from rmgpy.data.base import Entry, ForbiddenStructureException
 from rmgpy.data.thermo import *
 from rmgpy.data.solvation import *
 from rmgpy.data.kinetics import *
@@ -1282,10 +1282,25 @@ def addSeedMechanismToCore(self, seedMechanism, react=False):
         for entry in seedMechanism.entries.values():
             rxn = LibraryReaction(reactants=entry.item.reactants[:], products=entry.item.products[:], library=seedMechanism, kinetics=entry.data)
             r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
+            
+        # Perform species constraints and forbidden species checks
+        
         for spec in self.newSpeciesList:
+            if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
+                if 'allowed' in self.speciesConstraints and 'seed mechanisms' in self.speciesConstraints['allowed']:
+                    logging.warning("Species {0} from seed mechanism {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
+                else:
+                    raise ForbiddenStructureException("Species {0} from seed mechanism {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
+            if self.failsSpeciesConstraints(spec):
+                if 'allowed' in self.speciesConstraints and 'seed mechanisms' in self.speciesConstraints['allowed']:
+                    self.speciesConstraints['explicitlyAllowedMolecules'].append(spec.molecule[0])
+                    pass
+                else:
+                    raise ForbiddenStructureException("Species constraints forbids species {0} from seed mechanism {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, seedMechanism.label))
+
+        for spec in self.newSpeciesList:            
             if spec.reactive: spec.generateThermoData(database, quantumMechanics=self.quantumMechanics)
             spec.generateTransportData(database)
-        for spec in self.newSpeciesList:
             self.addSpeciesToCore(spec)
 
         for rxn in self.newReactionList:
@@ -1332,10 +1347,24 @@ def addReactionLibraryToEdge(self, reactionLibrary):
             rxn = LibraryReaction(reactants=entry.item.reactants[:], products=entry.item.products[:], library=reactionLibrary, kinetics=entry.data)
             r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
             if not isNew: logging.info("This library reaction was not new: {0}".format(rxn))
+            
+        # Perform species constraints and forbidden species checks
+        for spec in self.newSpeciesList:
+            if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
+                if 'allowed' in self.speciesConstraints and 'reaction libraries' in self.speciesConstraints['allowed']:
+                    logging.warning("Species {0} from reaction library {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
+                else:
+                    raise ForbiddenStructureException("Species {0} from reaction library {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
+            if self.failsSpeciesConstraints(spec):
+                if 'allowed' in self.speciesConstraints and 'reaction libraries' in self.speciesConstraints['allowed']:
+                    self.speciesConstraints['explicitlyAllowedMolecules'].append(spec.molecule[0])
+                    pass
+                else:
+                    raise ForbiddenStructureException("Species constraints forbids species {0} from reaction library {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, reactionLibrary.label))
+       
         for spec in self.newSpeciesList:
             if spec.reactive: spec.generateThermoData(database, quantumMechanics=self.quantumMechanics)
             spec.generateTransportData(database)
-        for spec in self.newSpeciesList:
             self.addSpeciesToEdge(spec)
 
         for rxn in self.newReactionList:
@@ -1676,3 +1705,41 @@ def saveChemkinFile(self, path, verbose_path, dictionaryPath=None, transportPath
             saveChemkinFile(verbose_path, speciesList, rxnList, verbose = True, checkForDuplicates=False)
             if dictionaryPath:
                 saveSpeciesDictionary(dictionaryPath, speciesList)
+                
+    def failsSpeciesConstraints(self, species):
+        """
+        Checks whether the species passes the speciesConstraints set by the user.  If not,
+        returns `True` for failing speciesConstraints.
+        """
+        explicitlyAllowedMolecules = self.speciesConstraints.get('explicitlyAllowedMolecules')
+        maxCarbonAtoms = self.speciesConstraints.get('maximumCarbonAtoms', 1000000)
+        maxHydrogenAtoms = self.speciesConstraints.get('maximumHydrogenAtoms', 1000000)
+        maxOxygenAtoms = self.speciesConstraints.get('maximumOxygenAtoms', 1000000)
+        maxNitrogenAtoms = self.speciesConstraints.get('maximumNitrogenAtoms', 1000000)
+        maxSiliconAtoms = self.speciesConstraints.get('maximumSiliconAtoms', 1000000)
+        maxSulfurAtoms = self.speciesConstraints.get('maximumSulfurAtoms', 1000000)
+        maxHeavyAtoms = self.speciesConstraints.get('maximumHeavyAtoms', 1000000)
+        maxRadicals = self.speciesConstraints.get('maximumRadicalElectrons', 1000000)
+        
+        struct = species.molecule[0]
+        for molecule in explicitlyAllowedMolecules:
+            if struct.isIsomorphic(molecule):
+                return False        
+        H = struct.getNumAtoms('H')
+        if struct.getNumAtoms('C') > maxCarbonAtoms:
+            return True
+        if H > maxHydrogenAtoms:
+            return True
+        if struct.getNumAtoms('O') > maxOxygenAtoms:
+            return True
+        if struct.getNumAtoms('N') > maxNitrogenAtoms:
+            return True
+        if struct.getNumAtoms('Si') > maxSiliconAtoms:
+            return True
+        if struct.getNumAtoms('S') > maxSulfurAtoms:
+            return True
+        if len(struct.atoms) - H > maxHeavyAtoms:
+            return True
+        if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H > 1):
+            return True
+        return False