trying to get rmg to get labled adj lists, and this is unfinished. se…

…e issue cfgoldsmith#72

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

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

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):

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)