diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 641a395db44..eb6905091be 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -32,6 +32,8 @@ This module contains functionality for working with kinetics families. """ +from __future__ import print_function +emily = open("label.txt", "w") # prints labels to a txt file in the top-level RMGpy folder import os.path import numpy as np import logging @@ -1220,6 +1222,12 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): # Also copy structures so we don't modify the originals # Since the tagging has already occurred, both the reactants and the # products will have tags + print("------------------", file=emily) + print("reactant structures", file=emily) + for reactant in reactantStructures: + print(reactant, file=emily) + print(reactant.toAdjacencyList(), file=emily) + if isinstance(reactantStructures[0], Group): reactantStructure = Group() else: @@ -1484,6 +1492,11 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): productStructures = [s for _, s in sorted(zip(lowest_labels, productStructures))] # Return the product structures + + print(label, file=emily) + print(reactantStructure.toAdjacencyList(label="reactant structure"), file=emily) + if reactantStructure.toAdjacencyList() != productStructure.toAdjacencyList(): + print(productStructure.toAdjacencyList(label="product structure"), file=emily) return productStructures def __generateProductStructures(self, reactantStructures, maps, forward): @@ -2122,6 +2135,13 @@ def generate_products_and_reactions(order): if isomorphic_species_lists(products, products0): rxnList.append(reaction) + # get reaction, if applicable + try: + print(rxnList[0], file=emily) + except: + pass + + # Determine the reactant-product pairs to use for flux analysis # Also store the reaction template (useful so we can easily get the kinetics later) for reaction in rxnList: @@ -2139,9 +2159,9 @@ def generate_products_and_reactions(order): # Unlabel the atoms for both reactants and products for species in itertools.chain(reaction.reactants, reaction.products): species.clearLabeledAtoms() - + # We're done with the labeled atoms, so delete the attribute - del reaction.labeledAtoms + # del reaction.labeledAtoms # Mark reaction reversibility reaction.reversible = self.reversible diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 2ccb81e4899..22fa408a944 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -605,7 +605,7 @@ def pressureDependence( rmg.pressureDependence.rmgmode = True def options(name='Seed', generateSeedEachIteration=False, saveSeedToDatabase=False, units='si', saveRestartPeriod=None, - generateOutputHTML=False, generatePlots=False, saveSimulationProfiles=False, verboseComments=False, + generateOutputHTML=False, generateLabeledAdjList=False, generatePlots=False, saveSimulationProfiles=False, verboseComments=False, saveEdgeSpecies=False, keepIrreversible=False, trimolecularProductReversible=True, wallTime='00:00:00:00'): rmg.name = name rmg.generateSeedEachIteration=generateSeedEachIteration @@ -614,7 +614,10 @@ def options(name='Seed', generateSeedEachIteration=False, saveSeedToDatabase=Fal rmg.saveRestartPeriod = Quantity(saveRestartPeriod) if saveRestartPeriod else None if generateOutputHTML: logging.warning('Generate Output HTML option was turned on. Note that this will slow down model generation.') - rmg.generateOutputHTML = generateOutputHTML + if generateLabeledAdjList: + logging.warning('Generate Labeled Adj List option was turned on.') + rmg.generateOutputHTML = generateOutputHTML + rmg.generateLabeledAdjList = generateLabeledAdjList rmg.generatePlots = generatePlots rmg.saveSimulationProfiles = saveSimulationProfiles rmg.verboseComments = verboseComments diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 2bf5de064cf..6d445eae6fe 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -70,7 +70,7 @@ import rmgpy.util as util from rmgpy.chemkin import ChemkinWriter -from rmgpy.rmg.output import OutputHTMLWriter +from rmgpy.rmg.output import OutputHTMLWriter, OutputLabeledAdjListWriter from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter from rmgpy.restart import RestartWriter from rmgpy.qm.main import QMDatabaseWriter @@ -123,6 +123,7 @@ class RMG(util.Subject): `saveRestartPeriod` The time period to periodically save a restart file (:class:`Quantity`), or ``None`` for never. `units` The unit system to use to save output files (currently must be 'si') `generateOutputHTML` ``True`` to draw pictures of the species and reactions, saving a visualized model in an output HTML file. ``False`` otherwise + `generateLabeledAdjList` ``True`` to export labeled adjacency lists in a text file `generatePlots` ``True`` to generate plots of the job execution statistics after each iteration, ``False`` otherwise `verboseComments` ``True`` to keep the verbose comments for database estimates, ``False`` otherwise `saveEdgeSpecies` ``True`` to save chemkin and HTML files of the edge species, ``False`` otherwise @@ -194,6 +195,7 @@ def clear(self): self.saveRestartPeriod = None self.units = 'si' self.generateOutputHTML = None + self.generateLabeledAdjList = None self.generatePlots = None self.saveSimulationProfiles = None self.verboseComments = None @@ -578,6 +580,9 @@ def register_listeners(self): if self.generateOutputHTML: self.attach(OutputHTMLWriter(self.outputDirectory)) + if self.generateLabeledAdjList: + self.attach(OutputLabeledAdjListWriter(self.outputDirectory)) + if self.saveRestartPeriod: warnings.warn("The option saveRestartPeriod is no longer supported and may be" " removed in version 2.3.", DeprecationWarning) @@ -594,7 +599,7 @@ def register_listeners(self): reactionSystem.attach(SimulationProfileWriter( self.outputDirectory, index, self.reactionModel.core.species)) reactionSystem.attach(SimulationProfilePlotter( - self.outputDirectory, index, self.reactionModel.core.species)) + self.outputDirectory, index, self.reactionModel.core.species)) def execute(self, **kwargs): diff --git a/rmgpy/rmg/output.py b/rmgpy/rmg/output.py index 75cbc3b0342..4a626a0e2ee 100644 --- a/rmgpy/rmg/output.py +++ b/rmgpy/rmg/output.py @@ -86,6 +86,7 @@ def saveOutputHTML(path, reactionModel, partCoreEdge='core'): if match: spec.index = int(match.group(0)[1:-1]) spec.label = spec.label[0:match.start()] + # Draw molecules if necessary fstr = os.path.join(dirname, 'species', '{0}.png'.format(spec)) if not os.path.exists(fstr): @@ -1314,7 +1315,7 @@ def saveOutput(rmg): """ logging.info('Saving current model core to HTML file...') saveOutputHTML(os.path.join(rmg.outputDirectory, 'output.html'), rmg.reactionModel, 'core') - + if rmg.saveEdgeSpecies == True: logging.info('Saving current model edge to HTML file...') saveOutputHTML(os.path.join(rmg.outputDirectory, 'output_edge.html'), rmg.reactionModel, 'edge') @@ -1346,3 +1347,185 @@ def __init__(self, outputDirectory=''): def update(self, rmg): saveOutput(rmg) + +def saveLabeledAdjList(path, reactionModel, partCoreEdge='core'): + + from rmgpy.rmg.model import PDepReaction + + if partCoreEdge == 'core': + reactions = [rxn for rxn in reactionModel.core.reactions] + reactionModel.outputReactionList + elif partCoreEdge == 'edge': + reactions = [rxn for rxn in reactionModel.edge.reactions] + reactionModel.outputReactionList + + familyCount = {} + for rxn in reactions: + + if isinstance(rxn, PDepReaction): + family = "PDepNetwork" + else: + family = rxn.getSource() + if family in familyCount: + familyCount[family] += 1 + else: + familyCount[family] = 1 + + + families_of_interest = familyCount.keys() + families_of_interest.sort() + import rmgpy + from rmgpy.molecule import Molecule + from rmgpy.species import Species + from rmgpy.reaction import Reaction + from rmgpy.data.rmg import RMGDatabase + rmg_database = RMGDatabase() + database_path = rmgpy.settings['database.directory'] + + # todo: make this work for all kinetic libraries + # this will not work for any other kinetic libraries that are multi level + libraries = os.listdir(os.path.join(database_path, "kinetics", "libraries" )) + cpox_libraries = ['Surface/CPOX_Pt/Deutschmann2006'] + + families_of_interest = [family for family in families_of_interest if family not in libraries] + families_of_interest = [family for family in families_of_interest if family not in cpox_libraries] + libraries_of_interest = [family for family in families_of_interest if family in libraries] + print libraries_of_interest + + rmg_database.load( + database_path, + kineticsFamilies=families_of_interest, + transportLibraries=[], + reactionLibraries=libraries_of_interest, + seedMechanisms=[], + thermoLibraries=[ + 'primaryThermoLibrary', + 'thermo_DFT_CCSDTF12_BAC', + 'CBS_QB3_1dHR'], + ) + + families = rmg_database.kinetics.families + write_string = "" + for reaction in reactions: + match = False + rmg_reactants = [] + rmg_products = [] + for react in reaction.reactants: + if isinstance(react, Species): + rmg_reactants.append(react.molecule) + elif isinstance(react, Molecule): + rmg_reactants.append([react]) + + for prod in reaction.products: + if isinstance(prod, Species): + rmg_products.append(prod.molecule) + elif isinstance(prod, Molecule): + rmg_products.append([prod]) + + test_reactants = [] + test_products = [] + if len(rmg_reactants) == 1: + test_reactants = rmg_reactants + elif len(rmg_reactants) == 2: + l1, l2 = rmg_reactants + for i in l1: + for j in l2: + test_reactants.append([i, j]) + elif len(rmg_products) == 3: + l1, l2, l3 = rmg_reactants + for i in l1: + for j in l2: + for k in l3: + test_reactants.append([i, j, k]) + + if len(rmg_products) == 1: + test_reactants = rmg_products + elif len(rmg_products) == 2: + l1, l2 = rmg_products + for i in l1: + for j in l2: + test_products.append([i, j]) + elif len(rmg_products) == 3: + l1, l2, l3 = rmg_products + for i in l1: + for j in l2: + for k in l3: + test_products.append([i, j, k]) + for name, family in families.items(): + if match: + continue + + labeled_reactants, labeled_products = (None, None) + for test_reactant in test_reactants: + for test_products in test_products: + if (labeled_reactants and labeled_products): + continue + labeled_reactants, labeled_products = family.getLabeledReactantsAndProducts( + reactants=test_reactant, + products=test_products + ) + + if (labeled_reactants) and (labeled_products): + print "{} matches {} reaction family.".format(reaction.__repr__(), name) + print + reactant_complex = Molecule() + write_string += reaction.__repr__() + write_string += "\n" + write_string += name #this is the string of the family name + write_string += "\n" + for reactant in labeled_reactants: + reactant.updateMultiplicity() + write_string += reactant.__repr__() + write_string += "\n" + write_string += reactant.toAdjacencyList() + write_string += "\n" + + reactant_complex = reactant_complex.merge(reactant) + reactant_complex.updateMultiplicity() + write_string += reactant_complex.__repr__() + write_string += "\n" + write_string += reactant_complex.toAdjacencyList() + write_string += "\n" + + product_complex = Molecule() + for product in labeled_products: + product.updateMultiplicity() + write_string += product.__repr__() + write_string += "\n" + write_string += product.toAdjacencyList() + write_string += "\n" + product_complex = product_complex.merge(product) + product_complex.updateMultiplicity() + write_string += product_complex.__repr__() + write_string += "\n" + write_string += product_complex.toAdjacencyList() + write_string += "\n" + match = True + if not match: + print "WARNING: Could not match {} to any provided reaction family...".format(reaction.__repr__()) + print + + f = open(os.path.join(path, "adjacency_lists.txt"), "w") + f.write(write_string) + f.close() + +def saveOutputAdj(rmg): + """ + Save the current reaction model to a pretty HTML file. + """ + logging.info('Saving current adj lists to file...') + saveLabeledAdjList(os.path.join(rmg.outputDirectory), rmg.reactionModel, 'core') + + +class OutputLabeledAdjListWriter(object): + """ + This class exports some labeled adj list + + Let's just copy as much of the write output html as possible + + Issue #72 in Franklin's rmg-cat + """ + def __init__(self, outputDirectory=''): + super(OutputLabeledAdjListWriter, self).__init__() + makeOutputSubdirectory(outputDirectory, 'species') + + def update(self, rmg): + saveOutputAdj(rmg) \ No newline at end of file