Merge pull request #1658 from ReactionMechanismGenerator/pdep_fix

Fix pressure dependent network generation in RMG jobs

arkane/explorerTest.py

@@ -76,7 +76,7 @@ def test_reactions(self):
         """
         test that the right number of reactions are in output network
         """
-        self.assertEqual(len(self.explorerjob.networks[0].pathReactions), 5)
+        self.assertEqual(len(self.explorerjob.networks[0].pathReactions), 6)
 
     def test_isomers(self):
         """

rmgpy/rmg/model.py

@@ -32,40 +32,30 @@
 Contains classes for working with the reaction model generated by RMG.
 """
 
-import logging
-import numpy
-import itertools
 import gc
+import itertools
+import logging
 import os
 
-import resource
-import psutil
-
-from sys import platform
-from multiprocessing import Pool
+import numpy
 
-from rmgpy.display import display
+import rmgpy.data.rmg
 from rmgpy import settings
 from rmgpy.constraints import failsSpeciesConstraints
-from rmgpy.quantity import Quantity
-from rmgpy.species import Species
-from rmgpy.molecule.molecule import Molecule
-from rmgpy.thermo.thermoengine import submit
-from rmgpy.reaction import Reaction
-from rmgpy.exceptions import ForbiddenStructureException
 from rmgpy.data.kinetics.depository import DepositoryReaction
 from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction
 from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction
-
-from rmgpy.kinetics import KineticsData, Arrhenius
-
 from rmgpy.data.rmg import getDB
-
-import rmgpy.data.rmg
-from .react import react_all
-from rmgpy.data.kinetics.common import ensure_independent_atom_ids, find_degenerate_reactions
+from rmgpy.display import display
+from rmgpy.exceptions import ForbiddenStructureException
+from rmgpy.kinetics import KineticsData, Arrhenius
+from rmgpy.quantity import Quantity
+from rmgpy.reaction import Reaction
+from rmgpy.rmg.pdep import PDepReaction, PDepNetwork
+from rmgpy.rmg.react import react_all
+from rmgpy.species import Species
+from rmgpy.thermo.thermoengine import submit
 
-from pdep import PDepReaction, PDepNetwork
 
 ################################################################################
 
@@ -599,14 +589,19 @@ def enlarge(self, newObject=None, reactEdge=False,
                 numOldEdgeReactions -= len(reactionsMovedFromEdge)
 
         else:
-            # We are reacting the edge
-            rxns = react_all(self.core.species, numOldCoreSpecies,
-                             unimolecularReact, bimolecularReact, trimolecularReact=trimolecularReact, procnum=procnum)
-
-            spcs = [self.retrieve_species(rxn) for rxn in rxns]
-
-            for rxn, spc in zip(rxns, spcs):
-                self.processNewReactions([rxn], spc, generateThermo=False)
+            # Generate reactions between all core species which have not been
+            # reacted yet and exceed the reaction filter thresholds
+            rxnLists, spcsTuples = react_all(self.core.species, numOldCoreSpecies,
+                                             unimolecularReact, bimolecularReact,
+                                             trimolecularReact=trimolecularReact,
+                                             procnum=procnum)
+
+            for rxnList, spcTuple in zip(rxnLists, spcsTuples):
+                if rxnList:
+                    # Identify a core species which was used to generate the reaction
+                    # This is only used to determine the reaction direction for processing
+                    spc = spcTuple[0]
+                    self.processNewReactions(rxnList, spc, generateThermo=False)
 
         ################################################################
         # Begin processing the new species and reactions
@@ -1865,28 +1860,6 @@ def retrieve(self, family_label, key1, key2):
         except KeyError: # no such short-list: must be new, unless in seed.
             return []
 
-    def getSpecies(self, obj):
-        """
-        Retrieve species object, by
-        polling the index species dictionary.
-        """
-        if isinstance(obj, int):
-            spc = self.indexSpeciesDict[obj]
-            return spc
-        return obj
-
-    def retrieve_species(self, rxn):
-        """
-        Searches for the first reactant or product in the reaction that is
-        a core species, which was used to generate the reaction in the first
-        place. Reactants or products not represented in the core will be
-        a newly-generated structure.
-        """
-        for obj in itertools.chain(rxn.reactants, rxn.products):
-            for spc in self.core.species:
-                if obj.isIsomorphic(spc):
-                    return spc
-        raise Exception("No core species were found in either reactants or products of {0}!".format(rxn))
 
 def generateReactionKey(rxn, useProducts=False):
     """

rmgpy/rmg/modelTest.py

@@ -28,9 +28,13 @@
 #                                                                             #
 ###############################################################################
 
+import itertools
 import os
 import unittest 
 
+import numpy as np
+from nose.plugins.attrib import attr
+
 from rmgpy import settings
 from rmgpy.data.rmg import RMGDatabase, database
 from rmgpy.rmg.main import RMG
@@ -118,7 +122,7 @@ def setUpClass(cls):
         for family in rmg.database.kinetics.families.values():
             family.forbidden = ForbiddenStructures()
         rmg.database.forbiddenStructures = ForbiddenStructures()
-    
+
     def testAddNewSurfaceObjects(self):
         """
         basic test that surface movement object management works properly
@@ -142,8 +146,8 @@ class item:
         spcs = [Species().fromSMILES('CC'), Species().fromSMILES('[CH3]')]
         spcTuples = [((spcA, spc), ['H_Abstraction']) for spc in spcs]
 
-        rxns = list(react(spcTuples, procnum))
-        rxns += list(react([((spcs[0], spcs[1]), ['H_Abstraction'])], procnum))
+        rxns = list(itertools.chain.from_iterable(react(spcTuples, procnum)))
+        rxns += list(itertools.chain.from_iterable(react([((spcs[0], spcs[1]), ['H_Abstraction'])], procnum)))
 
         for rxn in rxns:
             cerm.makeNewReaction(rxn)
@@ -246,7 +250,7 @@ def testMakeNewReaction(self):
         spcs = [Species().fromSMILES('CC'), Species().fromSMILES('[CH3]')]
         spcTuples = [((spcA, spc), ['H_Abstraction']) for spc in spcs]
 
-        rxns = list(react(spcTuples, procnum))
+        rxns = list(itertools.chain.from_iterable(react(spcTuples, procnum)))
 
         cerm = CoreEdgeReactionModel()
 
@@ -578,5 +582,123 @@ def tearDownClass(cls):
         rmgpy.data.rmg.database = None
 
 
+@attr('functional')
+class TestEnlarge(unittest.TestCase):
+    """
+    Contains unit tests for CoreEdgeReactionModel.enlarge.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        """
+        A method that is run ONCE before all unit tests in this class.
+        """
+        cls.dirname = os.path.abspath(os.path.join(os.path.dirname(__file__), 'temp'))
+        print cls.dirname
+        os.makedirs(os.path.join(cls.dirname, 'pdep'))
+
+        TESTFAMILY = 'R_Recombination'
+
+        cls.rmg = RMG()
+
+        from rmgpy.rmg.input import setGlobalRMG, pressureDependence
+        setGlobalRMG(cls.rmg)
+
+        pressureDependence(
+            method='modified strong collision',
+            maximumGrainSize=(0.5, 'kcal/mol'),
+            minimumNumberOfGrains=250,
+            temperatures=(300, 2100, 'K', 8),
+            pressures=(0.1, 100, 'bar', 5),
+            interpolation=('Chebyshev', 6, 4),
+            maximumAtoms=10,
+        )
+
+        cls.rmg.outputDirectory = cls.rmg.pressureDependence.outputFile = cls.dirname
+
+        cls.rmg.database = RMGDatabase()
+        cls.rmg.database.load(
+            path=settings['database.directory'],
+            thermoLibraries=['primaryThermoLibrary'],
+            kineticsFamilies=[TESTFAMILY],
+            reactionLibraries=[],
+        )
+
+        cls.rmg.reactionModel = CoreEdgeReactionModel()
+        cls.rmg.reactionModel.pressureDependence = cls.rmg.pressureDependence
+
+    def test_enlarge_1_add_nonreactive_species(self):
+        """Test that we can add a nonreactive species to CERM"""
+        m0 = Molecule(SMILES='[He]')
+        spc0 = self.rmg.reactionModel.makeNewSpecies(m0, label='He', reactive=False)[0]
+        self.rmg.reactionModel.enlarge(spc0)
+
+        self.assertEqual(len(self.rmg.reactionModel.core.species), 1)
+        self.assertFalse(self.rmg.reactionModel.core.species[0].reactive)
+
+    def test_enlarge_2_add_reactive_species(self):
+        """Test that we can add reactive species to CERM"""
+        m1 = Molecule(SMILES='CC')
+        spc1 = self.rmg.reactionModel.makeNewSpecies(m1, label='C2H4')[0]
+        self.rmg.reactionModel.enlarge(spc1)
+
+        self.assertEqual(len(self.rmg.reactionModel.core.species), 2)
+        self.assertTrue(self.rmg.reactionModel.core.species[1].reactive)
+
+        m2 = Molecule(SMILES='[CH3]')
+        spc2 = self.rmg.reactionModel.makeNewSpecies(m2, label='CH3')[0]
+        self.rmg.reactionModel.enlarge(spc2)
+
+        self.assertEqual(len(self.rmg.reactionModel.core.species), 3)
+        self.assertTrue(self.rmg.reactionModel.core.species[2].reactive)
+
+    def test_enlarge_3_react_edge(self):
+        """Test that enlarge properly generated reactions"""
+        self.rmg.reactionModel.enlarge(
+            reactEdge=True,
+            unimolecularReact=np.array([0, 1, 0], bool),
+            bimolecularReact=np.zeros((3, 3), bool),
+        )
+
+        self.assertEqual(len(self.rmg.reactionModel.edge.species), 2)
+        smiles = set([spc.SMILES for spc in self.rmg.reactionModel.edge.species])
+        self.assertEqual(smiles, {'[H]', 'C[CH2]'})
+
+        # We expect either C-C bond scission to be in the core and C-H bond scission to be in the edge
+        self.assertEqual(len(self.rmg.reactionModel.core.reactions), 1)
+        rxn = self.rmg.reactionModel.core.reactions[0]
+        smiles = set([spc.SMILES for spc in rxn.reactants + rxn.products])
+        self.assertEqual(smiles, {'CC', '[CH3]'})
+
+        self.assertEqual(len(self.rmg.reactionModel.edge.reactions), 1)
+        rxn = self.rmg.reactionModel.edge.reactions[0]
+        smiles = set([spc.SMILES for spc in rxn.reactants + rxn.products])
+        self.assertEqual(smiles, {'CC', 'C[CH2]', '[H]'})
+
+    def test_enlarge_4_create_pdep_network(self):
+        """Test that enlarge properly creates a pdep network"""
+        self.assertEqual(len(self.rmg.reactionModel.networkList), 1)
+        self.assertEqual(len(self.rmg.reactionModel.networkList[0].source), 1)
+        self.assertEqual(self.rmg.reactionModel.networkList[0].source[0].label, 'C2H4')
+
+        self.assertEqual(len(self.rmg.reactionModel.networkDict), 1)
+        self.assertEqual(len(self.rmg.reactionModel.networkDict.keys()[0]), 1)
+        self.assertEqual(self.rmg.reactionModel.networkDict.keys()[0][0].label, 'C2H4')
+
+    @classmethod
+    def tearDownClass(cls):
+        """
+        A method that is run ONCE after all unit tests in this class.
+
+        Clear global variables and clean up files.
+        """
+        import rmgpy.data.rmg
+        rmgpy.data.rmg.database = None
+        import rmgpy.rmg.input
+        rmgpy.rmg.input.rmg = None
+        import shutil
+        shutil.rmtree(cls.dirname)
+
+
 if __name__ == '__main__':
     unittest.main()

rmgpy/rmg/parreactTest.py

@@ -32,6 +32,7 @@
 This module contains unit tests of the rmgpy.parallel module.
 """
 
+import itertools
 import os
 import sys
 import unittest
@@ -93,7 +94,7 @@ def generate():
     spcTuples = [(spcA, spc) for spc in spcs]
     procnum = 2
 
-    reactionList = list(react(spcTuples, procnum))
+    reactionList = list(itertools.chain.from_iterable(react(spcTuples, procnum)))
 
     if not reactionList: return False