diff --git a/arkane/explorerTest.py b/arkane/explorerTest.py index f95bd15309..9245f49b6c 100644 --- a/arkane/explorerTest.py +++ b/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): """ diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 8187362050..18be219271 100644 --- a/rmgpy/rmg/model.py +++ b/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): """ diff --git a/rmgpy/rmg/modelTest.py b/rmgpy/rmg/modelTest.py index b81c606584..097a320726 100644 --- a/rmgpy/rmg/modelTest.py +++ b/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() diff --git a/rmgpy/rmg/parreactTest.py b/rmgpy/rmg/parreactTest.py index a6b360ca41..74b8fad9ad 100644 --- a/rmgpy/rmg/parreactTest.py +++ b/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 diff --git a/rmgpy/rmg/react.py b/rmgpy/rmg/react.py index db901702ca..48dc251161 100644 --- a/rmgpy/rmg/react.py +++ b/rmgpy/rmg/react.py @@ -31,43 +31,47 @@ """ Contains functions for generating reactions. """ -import itertools import logging from rmgpy.data.rmg import getDB from multiprocessing import Pool ################################################################################ -def react(spc_tuples, procnum=1): + + +def react(spc_fam_tuples, procnum=1): """ - Generate reactions between the species in the - list of species tuples for all the reaction families available. + Generate reactions between the species in the list of species-family tuples + for the optionally specified reaction families. - For each tuple of one or more Species objects [(spc1,), (spc2, spc3), ...] - the following is done: + Each item should be a tuple of a species tuple and an optional family list: + [((spc1,), [family1, ...]), ((spc2, spc3), [family2, ...]), ...] + OR + [((spc1,),), ((spc2, spc3),), ...] - A list of tuples is created for each resonance isomer of the species. - Each tuple consists of (Molecule, index) with the index the species index of the Species object. + If no family list is provided, all of the loaded families are considered. - Possible combinations between the first spc in the tuple, and the second species in the tuple - is obtained by taking the combinatorial product of the two generated [(Molecule, index)] lists. + Args: + spc_fam_tuples (list): list of tuples for reaction generation + procnum (int, optional): number of processors used for reaction generation - Returns a flat generator object containing the generated Reaction objects. + Returns: + list of lists of reactions generated from each species tuple (note: empty lists are possible) """ # Execute multiprocessing map. It blocks until the result is ready. # This method chops the iterable into a number of chunks which it # submits to the process pool as separate tasks. if procnum == 1: logging.info('For reaction generation {0} process is used.'.format(procnum)) - reactions = map(_react_species_star, spc_tuples) + reactions = map(_react_species_star, spc_fam_tuples) else: logging.info('For reaction generation {0} processes are used.'.format(procnum)) p = Pool(processes=procnum) - reactions = p.map(_react_species_star, spc_tuples) + reactions = p.map(_react_species_star, spc_fam_tuples) p.close() p.join() - return itertools.chain.from_iterable(reactions) + return reactions def _react_species_star(args): @@ -79,6 +83,13 @@ def react_species(species_tuple, only_families=None): """ Given a tuple of Species objects, generates all possible reactions from the loaded reaction families and combines degenerate reactions. + + Args: + species_tuple (tuple): tuple of 1-3 Species objects to react together + only_families (list, optional): list of reaction families to consider + + Returns: + list of generated reactions """ species_tuple = tuple([spc.copy(deep=True) for spc in species_tuple]) @@ -90,8 +101,22 @@ def react_species(species_tuple, only_families=None): def react_all(core_spc_list, numOldCoreSpecies, unimolecularReact, bimolecularReact, trimolecularReact=None, procnum=1): """ - Reacts the core species list via uni-, bi-, and trimolecular - reactions and splits reaction families per task for improved load balancing in parallel runs. + Reacts the core species list via uni-, bi-, and trimolecular reactions. + + For parallel processing, reaction families are split per task for improved + load balancing. This is currently hard-coded using reaction family labels. + + Args: + core_spc_list (list): list of all core species + numOldCoreSpecies (int): current number of core species in the model + unimolecularReact (np.ndarray): reaction filter flags indicating which species to react unimolecularly + bimolecularReact (np.ndarray): reaction filter flags indicating which species to react bimolecularly + trimolecularReact (np.ndarray, optional): reaction filter flags indicating which species to react trimolecularly + procnum (int, optional): number of processors used for reaction generation + + Returns: + a list of lists of reactions generated from each species tuple + a list of species tuples corresponding to each list of reactions """ # Select reactive species that can undergo unimolecular reactions: spc_tuples = [(core_spc_list[i],) @@ -147,5 +172,5 @@ def react_all(core_spc_list, numOldCoreSpecies, unimolecularReact, bimolecularRe else: spc_fam_tuples.append((spc_tuple, )) - return list(react(spc_fam_tuples, procnum)) + return react(spc_fam_tuples, procnum), [fam_tuple[0] for fam_tuple in spc_fam_tuples] diff --git a/rmgpy/rmg/reactTest.py b/rmgpy/rmg/reactTest.py index e64c4452db..ab9fa8c57a 100644 --- a/rmgpy/rmg/reactTest.py +++ b/rmgpy/rmg/reactTest.py @@ -28,6 +28,7 @@ # # ############################################################################### +import itertools import os import unittest import numpy as np @@ -76,7 +77,7 @@ def testReact(self): spcs = [Species().fromSMILES('CC'), Species().fromSMILES('[CH3]')] spc_tuples = [((spc_a, spc), ['H_Abstraction']) for spc in spcs] - reaction_list = list(react(spc_tuples, procnum)) + reaction_list = list(itertools.chain.from_iterable(react(spc_tuples, procnum))) self.assertIsNotNone(reaction_list) self.assertEqual(len(reaction_list), 3) self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in reaction_list])) @@ -93,7 +94,7 @@ def testReactParallel(self): spcs = [Species().fromSMILES('CC'), Species().fromSMILES('[CH3]')] spc_tuples = [((spc_a, spc), ['H_Abstraction']) for spc in spcs] - reaction_list = list(react(spc_tuples, procnum)) + reaction_list = list(itertools.chain.from_iterable(react(spc_tuples, procnum))) self.assertIsNotNone(reaction_list) self.assertEqual(len(reaction_list), 3) self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in reaction_list])) @@ -112,10 +113,14 @@ def testReactAll(self): ] n = len(spcs) - reaction_list = react_all(spcs, n, np.ones(n), np.ones([n, n]), np.ones([n, n, n]), procnum) + reaction_list, spc_tuples = react_all(spcs, n, np.ones(n), np.ones([n, n]), np.ones([n, n, n]), procnum) self.assertIsNotNone(reaction_list) - self.assertEqual(len(reaction_list), 44) - self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in reaction_list])) + self.assertEqual(len(reaction_list), 34) + self.assertEqual(len(spc_tuples), 34) + + flat_rxn_list = list(itertools.chain.from_iterable(reaction_list)) + self.assertEqual(len(flat_rxn_list), 44) + self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in flat_rxn_list])) def testReactAllParallel(self): """ @@ -133,10 +138,14 @@ def testReactAllParallel(self): ] n = len(spcs) - reaction_list = react_all(spcs, n, np.ones(n), np.ones([n, n]), np.ones([n, n, n]), procnum) + reaction_list, spc_tuples = react_all(spcs, n, np.ones(n), np.ones([n, n]), np.ones([n, n, n]), procnum) self.assertIsNotNone(reaction_list) - self.assertEqual(len(reaction_list), 44) - self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in reaction_list])) + self.assertEqual(len(reaction_list), 94) + self.assertEqual(len(spc_tuples), 94) + + flat_rxn_list = list(itertools.chain.from_iterable(reaction_list)) + self.assertEqual(len(flat_rxn_list), 44) + self.assertTrue(all([isinstance(rxn, TemplateReaction) for rxn in flat_rxn_list])) def tearDown(self): """