Unified MEASURE input and output file formats.

As of now the only difference between the two is the presence of
pdepreaction() blocks in the latter, signifying the net reactions and their
pressure-dependent kinetics. This means that the entirety of the input file
is contained in the output file. Accordingly, there are no longer separate
readInput() and readOutput() methods, nor separate writeInput() and
writeOutput() methods; instead, there are now simply readFile() and
writeFile() methods to use for both cases.

Author: jwallen

rmgpy/measure/input.py

@@ -42,7 +42,7 @@
 from rmgpy.reaction import Reaction
 from rmgpy.molecule import Molecule
 from rmgpy.statmech import *
-from rmgpy.kinetics import Arrhenius
+from rmgpy.kinetics import Arrhenius, PDepArrhenius, Chebyshev
 from rmgpy.thermo import *
 
 from network import Network
@@ -141,6 +141,20 @@ def reaction(reactants, products, kinetics=None, transitionState=None):
     network.pathReactions.append(rxn)
     logging.debug('Found reaction "{0}"'.format(rxn))
 
+def pdepreaction(reactants, products, kinetics=None):
+    global network
+    try:
+        rxn = Reaction(
+            reactants = reactants,
+            products = products,
+            kinetics=kinetics,
+        )
+    except KeyError, e:
+        raise InputError('A reaction was encountered with species "{0}", but that species was not found in the input file.'.format(e.args[0]))
+        
+    network.netReactions.append(rxn)
+    logging.debug('Found pdepreaction "{0}"'.format(rxn))
+
 def collisionModel(type, parameters, bathGas):
     global network
     if type.lower() == 'single exponential down':
@@ -309,9 +323,9 @@ def getPressuresForModel(model, Pmin, Pmax, Pcount):
         Plist = 10.0 ** numpy.linspace(math.log10(Pmin), math.log10(Pmax), Pcount)
     return Plist
 
-def readInput(path):
+def readFile(path):
     """
-    Reads a MEASURE input file from location `path` on disk. The input file
+    Reads a MEASURE input or output file from location `path` on disk. The file
     format is described in the :ref:`measureusersguide`. Returns a number of
     quantities:
 
@@ -381,6 +395,9 @@ def readInput(path):
         'reaction': reaction,
         'Arrhenius': Arrhenius,
         'TransitionState': TransitionState,
+        'pdepreaction': pdepreaction,
+        'Chebyshev': Chebyshev,
+        'PDepArrhenius': PDepArrhenius,
         'collisionModel': collisionModel,
         'temperatures': temperatures,
         'pressures': pressures,
@@ -451,6 +468,11 @@ def readInput(path):
             rxn.reactants.sort()
             rxn.products = [speciesDict[label] for label in rxn.products]
             rxn.products.sort()
+        for rxn in network.netReactions:
+            rxn.reactants = [speciesDict[label] for label in rxn.reactants]
+            rxn.reactants.sort()
+            rxn.products = [speciesDict[label] for label in rxn.products]
+            rxn.products.sort()
 
         # Figure out which configurations are isomers, reactant channels, and product channels
         for rxn in network.pathReactions:

rmgpy/measure/main.py

@@ -127,8 +127,8 @@ def execute(inputFile, outputFile=None, drawFile=None, logFile=None, quiet=False
     logHeader()
 
     # Load input file
-    from rmgpy.measure.input import readInput
-    params = readInput(os.path.relpath(inputFile))
+    from rmgpy.measure.input import readFile
+    params = readFile(os.path.relpath(inputFile))
 
     # Only proceed if the input network is valid
     if params is not None:
@@ -188,12 +188,12 @@ def execute(inputFile, outputFile=None, drawFile=None, logFile=None, quiet=False
             logging.info('')
 
             # Save results to file
-            from rmgpy.measure.output import writeOutput
+            from rmgpy.measure.output import writeFile
             if outputFile is not None:
                 out = os.path.relpath(outputFile)
             else:
                 out = os.path.join(outputDirectory, 'output.py')
-            writeOutput(out, network, Tlist, Plist, Elist, method, model)
+            writeFile(out, network, Tlist, Plist, Elist, method, model, Tmin, Tmax, Pmin, Pmax)
 
     # Log end timestamp
     logging.info('')

rmgpy/measure/output.py

@@ -173,7 +173,6 @@ def writeReaction(f, rxn):
     f.write('reaction(\n')
     f.write('    reactants=[{0}],\n'.format(', '.join([('"{0}"'.format(spec)) for spec in rxn.reactants])))
     f.write('    products=[{0}],\n'.format(', '.join([('"{0}"'.format(spec)) for spec in rxn.products])))
-    f.write('    reversible={0!r},\n'.format(rxn.reversible))
     if rxn.kinetics is not None:
         f.write('    kinetics={0!r},\n'.format(rxn.kinetics))
     if rxn.transitionState is not None:
@@ -205,7 +204,6 @@ def writePDepReaction(f, rxn):
     f.write('pdepreaction(\n')
     f.write('    reactants=[{0}],\n'.format(', '.join([('"{0}"'.format(spec)) for spec in rxn.reactants])))
     f.write('    products=[{0}],\n'.format(', '.join([('"{0}"'.format(spec)) for spec in rxn.products])))
-    f.write('    reversible={0!r},\n'.format(rxn.reversible))
     if rxn.kinetics is not None:
         if isinstance(rxn.kinetics, Chebyshev):
             f.write('    kinetics=Chebyshev(\n')
@@ -232,67 +230,10 @@ def writePDepReaction(f, rxn):
 
 ################################################################################
 
-def writeOutput(path, network, Tlist, Plist, Elist, method, model):
-    """
-    Write a MEASURE output file to `path` on disk. The parameters needed mirror
-    those returned by :meth:`readInput()`:
-
-    * The :class:`Network` object `network` representing the unimolecular
-      reaction network
-
-    * The list of temperatures `Tlist` in K to be used in the master equation
-      calculation
-
-    * The list of pressures `Plist` in Pa to be used in the master equation
-      calculation
-
-    * A tuple `Elist` containing the maximum energy grain size in J/mol and the
-      minimum number of energy grains to use in the master equation calculation;
-      whichever of these results in more energy grains
-
-    * The approximate `method` to use to estimate the phenomenological rate
-      coefficients :math:`k(T,P)`
-
-    * The interpolation `model` to fit the estimated :math:`k(T,P)` values to
-
-    If successful, the file created on disk will contain all of the species
-    and net reaction data, including all phenomenological rate coefficients
-    :math:`k(T,P)`.
-    """
-
-    logging.info('Saving output to "{0}"...'.format(path))
-
-    f = open(os.path.relpath(path), 'w')
-
-    f.write('################################################################################\n')
-    f.write('#\n')
-    f.write('#   MEASURE output file for {0}\n'.format(network))
-    f.write('#\n')
-    f.write('#   Generated on {0}\n'.format(time.asctime()))
-    f.write('#\n')
-    f.write('################################################################################\n\n')
-    
-    # Write metadata
-    writeNetworkMetadata(f, network)
-    
-    # Write each species to file
-    writeNetworkSpecies(f, network)
-    f.write('################################################################################\n\n')
-    # Write each configuration to file
-    writeNetworkConfigurations(f, network)
-    f.write('################################################################################\n\n')
-    
-    # Write each net reaction to file
-    writeNetworkNetReactions(f, network)
-    
-    f.close()
-
-################################################################################
-
-def writeInput(path, network, Tlist, Plist, Elist, method, model):
+def writeFile(path, network, Tlist, Plist, Elist, method, model, Tmin, Tmax, Pmin, Pmax):
     """
     Write a MEASURE input file to `path` on disk. The parameters needed mirror
-    those returned by :meth:`readInput()`:
+    those returned by :meth:`readFile()`:
 
     * The :class:`Network` object `network` representing the unimolecular
       reaction network
@@ -320,7 +261,7 @@ def writeInput(path, network, Tlist, Plist, Elist, method, model):
 
     f.write('################################################################################\n')
     f.write('#\n')
-    f.write('#   MEASURE input file for {0}\n'.format(network))
+    f.write('#   MEASURE file for {0}\n'.format(network))
     f.write('#\n')
     f.write('#   Generated on {0}\n'.format(time.asctime()))
     f.write('#\n')
@@ -353,14 +294,26 @@ def writeInput(path, network, Tlist, Plist, Elist, method, model):
         f.write('    },\n')
     f.write(')\n\n')
 
-    f.write('temperatures(([{0}],"K"))\n'.format(', '.join([('{0:g}'.format(T)) for T in Tlist])))
-    f.write('pressures(([{0}],"bar"))\n'.format(', '.join([('{0:g}'.format(P/1e5)) for P in Plist])))
-    dE, count = Elist
+    if Tmin is not None and Tmax is not None:
+        f.write('temperatures(Tmin=({0},"K"), Tmax=({1},"K"), count={2})\n'.format(Tmin, Tmax, len(Tlist)))
+    else:
+        f.write('temperatures(([{0}],"K"))\n'.format(', '.join([('{0:g}'.format(T)) for T in Tlist])))
+    if Pmin is not None and Pmax is not None:
+        f.write('pressures(Pmin=({0},"bar"), Pmax=({1},"bar"), count={2})\n'.format(Pmin/1e5, Pmax/1e5, len(Plist)))
+    else:
+        f.write('pressures(([{0}],"bar"))\n'.format(', '.join([('{0:g}'.format(P/1e5)) for P in Plist])))
+    dE = Elist[1] - Elist[0]; count = len(Elist)
     f.write('energies(dE=({0:g},"kJ/mol"), count={1:d})\n'.format(dE/1000.0, count))
     f.write('method("{0}")\n'.format(method))
     if model[0].lower() == 'chebyshev':
         f.write('interpolationModel("chebyshev", {0:d}, {1:d})\n'.format(model[1], model[2]))
     else:
         f.write('interpolationModel("{0}")\n'.format(model[0].lower()))
-
+    f.write('\n')
+    
+    if len(network.netReactions) > 0:
+        f.write('################################################################################\n\n')
+        # Write each net reaction to file
+        writeNetworkNetReactions(f, network)
+    
     f.close()